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ABSTRACT 


A  new  method  for  the  numerical  simulation  of  three-dimensional 
incompressible  flows  is  described.  Our  vortex-in-cell  (VIC)  method 
traces  the  motion  of  the  vortex  filaments  in  the  velocity  field  which 
these  filaments  create.  The  velocity  field  is  not  calculated;  directly 
by  the  Biot-Savart  law  of  interaction  but^by  creating  a  mesh-record  of 
the  vorticity  field,  then  integrating  a  Poisson's  equation  via  the 
fast  Fourier  transform  to  get  the  stream  function  and  generating  a 
mesh-record  of  the  velocity  field.  The  computed  scales  of  motion  are 
assumed  to  be  essentially  inviscid.  Viscous  or  subgrid-scale  effects 
are  incorporated  into  a  filtering  procedure  in  wave  vector  space. 

Three  computational  experiments  were  pursued  in  three-dimensional 
space.  The  velocity  of  translation  of  a  single  vortex  ring  was 
measured  and  compared  with  the  Biot-Savart  law^j  The  agreement  was 
excellent. )»Next,  the  evolution  in  time  of  an  infinite  periodic  array 
of  closed  vortex  filaments  (Taylor-Green)  was  studied,  \  The  results 
were  compared  with  those  obtained  previously  by  the  spectral  method. 

The  Fourier  mode  spectra  were  in  very  good  agreement .V  The  third  simu¬ 
lation  follows  a  mixing  layer  from  an  initial  state  of  uniform  vorticity 
with  two-  and  three-dimensional  small  perturbations.  Streamwise  per¬ 
turbations  lead  to  the  usual  roll-up  of  vortex  patterns  with  spanwise 
uniformity.  Combined  with  spanwise  perturbations,  it  was  observed  that 

r\ _ 

the  third  dimension  is  extremely  important  in  the  evolution  of  a  mixing 
layer  even  though  it  is  an  "ignorable  coordinate"  in  the  idealized 
uniform  shear  flow.  Substantial  streamwise  vortex  deformation  and 
stretching  is  generated  by  the  spanwise  perturbations. 
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Chapter  I 
INTRODUCTION 

It  la  becoming  evident  that  a  proper  understanding  of  turbulence 
will  depend  heavily  on  numerical  methods  of  solving  the  Navier-Stokes 
equations  In  three  dimensions  at  high  Reynolds  numbers.  In  this  work, 
we  describe  a  new  method  for  the  numerical  simulation  of  three-dimen¬ 
sional  incompressible  flows.  Our  approach  differs  from  other  numerical 
fluid  simulations  in  that,  rather  than  solving  the  Navler-Stokes  equa¬ 
tions  on  an  Eulerian  mesh,  we  emphasize  the  vortical  part  of  the  flow 
by  solving  the  vorticity  equation,  using  a  hybrid  method. 

Numerics  means  discretization,  and  in  breaking  up  the  fluid  or  gas 
into  discrete  elements,  one  can,  and  should,  make  use  of  the  fact  that 
vortices  naturally  preserve  their  identity.  Even  when  the  vorticity 
is  itself  continuous,  as  in  laminar  shear  flow,  the  parametrization 
of  fluid  elements  according  to  (infinitesimal)  vortex  filaments  helps 
greatly  toward  programming  hydro-  or  gas-dynamics  into  the  computer. 

Historically,  the  first  numerical  calculation  using  a  two-dimen¬ 
sional  discrete  vortex  element  method  was  made  by  Rosenhead^1  ,  using 
a  few  vortices.  Since  then,  the  same  method  has  been  applied  to  various 
two-dimensional  shear  flows  and  is  summarized  in  a  literature  survey  by 

T2  T 

Fink  and  SohL  .  Most  of  these  efforts  were  directed  toward  understand¬ 
ing  the  time  evolution  of  finite-area  vorticity  regions  or  the  initial 
break-up  of  the  laminar  shear  layer  (usually  done  for  short  times  or 
small  regions  with  periodic  boundary  conditions).  Some  of  these  recent 
computer  calculations  involved  thousands  of  vortices^’^’^.  Other  impor¬ 
tant  flows  were  also  considered:  two-dimensional  turbulence^’^’®’^, 
separated  flows^0^,  and  stratified  flows  in  both  homogeneous  and  porous 


[11,12] 
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A  flow  field  can  be  represented  to  any  required  precision  by  a 
sufficiently  large  assembly  of  discrete  vortices  and  the  time  evolution 
of  the  field  can  be  transcribed  into  a  kind  of  particle  mechanics  of 
these  vortices.  If  one  integrates  analytically,  in  two  or  three  dimen¬ 
sions,  the  Poisson's  equation  relating  the  velocity  field  to  the  vor- 
ticity  field  (see  Eq.  2.4),  the  result  is  a  direct  Biot-Savart  integra- 

[34] 

tion  or  summation  over  all  the  vortices  making  the  vorticity  field 
With  N  vortices,  a  time  step  in  a  direct  summation  scheme  will  in¬ 
volve  evaluating  N-l  terms  for  the  velocity  which  displaces  a  single 
vortex,  and  thus  on  the  order  of  NXN  terms  per  time  step.  For  large 

N  ,  such  a  code  will  be  very  time  consuming.  Thus,  the  two-dimensional 

[7] 

mixing  layer  calculation  of  Ashurst  ,  where  the  number  of  vortices  was 
increased  from  one  to  eight  hundred,  required  two  hundred  and  fifty  hours 
of  computing  time  on  a  CDC  6600.  Understandably,  previous  direct  summa- 

[13] 

tion  calculations  were  much  more  modest.  Michalke  studied  the 
linear  and  partly  nonlinear  instability  of  seventy-two  vortices  arranged 

[i41 

on  three  close  parallel  lines.  Acton  ,  still  in  two  dimensions, 
studied  the  same  problem,  but  with  ninety-six  vortices  arranged  initially 
on  four  close  sinusoids,  and  ran  for  a  longer  time.  In  Acton's  simula¬ 
tion,  the  sinusoid  is  seen  to  amplify,  break  and  form  two  vortex  blobs 
which  then  merge  into  one. 

An  alternative  method  for  time-stepping  a  set  of  parallel  vortices 
in  two  dimensions,  the  cloud-in-cell  (CIC)  method  as  it  is  sometimes 
called,  was  described  by  Christiansen^-1"’’'’’1^ .  In  this  algorithm,  the 
basic  variables  are  still  the  positions  and  strengths  of  the  vortices, 
but  now  a  grid  is  laid  down  in  the  plane  perpendicular  to  the  vortex 
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filaments  and  covering  the  flow  area.  At  every  time  step,  a  grid  vor- 
ticity  is  generated  by  distributing  the  vorticity  from  each  vortex  over 
the  four  neighboring  grid  points  suitably  weighted.  This  grid  vorticity 
is  then  used  to  generate  a  grid  stream  function  by  solving  Poisson's 
equation.  The  stream  function  is  differenced  on  the  grid  to  produce  a 
velocity  field,  which  is  finally  interpolated  back  to  the  vortex  posi¬ 
tions.  The  advantage  of  this  apparently  rather  elaborate  detour  is 
that  the  Poisson  inversion  can  be  accomplished  by  fast  Fourier  transform 
techniques.  If  the  grid  is  MXM  ,  this  is  an  order  M^log^M  calcula¬ 
tion.  Since  in  a  typical  2-D  calculation  today,  the  number  of  grid 
squares  MXM  is  some  small  multiple  of  the  total  number  of  vortices  N  , 
this  required  on  the  order  of  N  log^N  operations  per  time  step.  The 
smoothing  and  interpolations  take  of  the  order  N  operations.  When  N 
is  large,  N  log^N  is  considerably  smaller  than  NXN  . 

A  disadvantage  of  the  cloud-in-cell  method  (as  against  the  direct 
interaction  method)  is  having  to  watch  out  for  grid  effects.  In  parti¬ 
cular,  "aliassing"  should  be  avoided  or  minimized,  as  in  the  case  of 

[28  57] 

conventional  Eulerian  fluid  simulations  .  To  counteract  grid 

r  171 

effects,  Wang  did  some  extensive  two-dimensional  simulations  in 

which  he  improved  on  the  CIC  method  by  first,  using  cubic  splines  for 

interpolation  (i.e.  referring  to  the  nearest  l6  grid  points  for  each 

vortex  node  rather  than  to  the  nearest  4)  and  second,  applying  a  gaus- 

slan  shape  factor  or  "filter"  in  wave  vector  space.  The  insensitivity 

of  the  resulting  potential  contours  or  flow  lines  to  the  underlying  grid 

Cl7] 

is  shown  in  some  diagrams  presented  inu  . 

The  "filter"  can  serve  a  dual  purpose  in  vortex  methods.  It 
obviously  helps  to  de-emphasize  the  high  harmonics  which  are  the  most 
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sensitive  to  grid  effects.  But  it  also  spreads  the  individual  vortex 
strands  from  singular  filaments  to  tubes  of  finite  thickness.  It  was 
pointed  out1  J  that,  because  high  induced  velocities  occur  when  two 
vortices  come  close  together,  unless  a  finite  core  radius  is  used,  the 
accuracy  of  the  discrete  vortex  element  method  does  not  improve  as  more 
and  more  elementary  vortices  are  employed  to  represent  a  given  vorti- 
city  field.  This  small  cut-off  radius  represents  the  smallest  scale 
we  can  compute.  The  effect  is  not  cumulative;  the  vorticity  is  spread 
a  little  way  only. 

It  is  well  known  in  turbulence  theory  that  energy  is  transferred 

from  large  to  small  scale  through  the  three-dimensional  vortex-stretching 

no! 

mechanism.  Although  many  explanations  have  been  given1-  7  ,  there  seems 
to  be  little  quantitative  support.  It  has  also  long  been  understood 
that  jets  and  wakes  can  be  idealized  by  a  series  of  vortex  rings,  but 
few  studies  have  been  carried  out  in  terms  of  three-dimensional  vorticity 
theory. 

In  three  dimensions,  the  direct  summation  scheme  has  been  applied 

[20  21  22  23] 

to  a  small  number  of  simple  vortex  filaments  »  »  »  j  but  at 

considerable  cost  because  of  the  time  required  to  sum  all  the  mutual 
Biot-Savart  interactions  between  the  many  elements  in  all  the  filaments. 
In  our  method,  the  velocity  field  is  not  calculated  directly  by  the  Biot- 
Savart  law  of  interaction  but  by  creating  a  mesh-record  of  the  vorticity 
field,  then  integrating  a  Poisson's  equation  to  generate  a  mesh-record 
of  the  velocity  field.  The  "cell"  or  "mesh"  method  speeds  up  the  cal¬ 
culations  of  the  interactions  and  allows  the  three-dimensional  vortex 
pushing  method  to  be  applied  to  a  space  densely  filled  with  vortex  fila¬ 
ments,  each  filament  being  resolved  in  fine  detail  along  its  length. 
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Some  features  of  our  method  are  the  following:  a)  each  vortex  of 
finite  size  has  a  gaussian  vorticity  distribution,  which  is  achieved  by 
combining  quadratic  spline  interpolation  with  spectral  shaping  in  three- 
dimensional  wave  vector  space  where  Poisson's  equation  for  the  velocity 
field  is  solved;  b)  a  machine-coded  FFT  is  applied  to  transform  vorti¬ 
city  records  from  real  space  to  wave  vector  space,  and  to  transform 
back  velocity  records  every  time  step;  c)  only  the  spline  coeffi¬ 
cients  of  the  physical  variables  are  accumulated  or  recorded  on  grid 
points.  Good  subgrid  accuracy  is  achieved  through  the  spline  inter¬ 
polation,  minimizing  grid  effects  and  preventing  numerical  instability. 
Some  of  the  basic  advantages  of  our  discrete  vortex  element  method  over 
the  finite-difference  methods  for  numerical  simulation  of  flows  are  as 
follows:  a)  the  method  lends  itself  to  simulation  of  flows  at  high 
Reynolds  number;  b)  it  eliminates  the  necessity  to  solve  for  the  irro- 
tational  flow  domain,  which  remains  passive;  c )  it  eliminates  the  Cour- 
ant-Friedrichs-Levy  (CFL)  condition,  which  imposes  a  time-step  limita¬ 
tion  based  upon  the  arbitrarily  chosen  mesh  size  instead  of  the  physical 
time  scale  based  upon  vortex  interactions. 

Basically,  three  computational  experiments  were  pursued.  The 
first  one  involves  vortex  rings.  The  velocity  of  translation  of  a 
single  vortex  ring  in  three-dimensional  space  was  measured  in  a  vortex- 
in-cell  (VIC)  run.  This  can  also  be  calculated  analytically,  at  some 
cost,  using  the  well-known  Green's  function  method  (a  Lagrangian  method 
using  the  Biot-Savart  law).  We  compare  the  results  of  the  two  methods, 
(see  Chapter  IV,  section  1). 

The  second  experiment  considers  a  particular  periodic  vortex  pat- 

r24l 

tern  due  to  Taylor  and  Greenu  .  After  attempting  a  perturbation 
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analysis  to  trace  the  evolution  of  this  flow,  it  appeared  that,  as  in 
Taylor  and  Green's  original  work,  such  an  analysis  would  fail  to  pro¬ 
vide  medium  to  large  time  solutions,  or  high  Reynolds  number  solutions. 
Our  vortex-in-cell  method  was  used  to  print  out  Fourier  mode  energy  and 
enstrophy  (mean  square  vorticity).  These  results  were  compared  with  a 
purely  spectral  (Galerkin)  method,  most  appropriate  to  simulate  such  a 
naturally  periodic  flow.  Also,  a  movie  tracing  the  vortex  lines  in 
three-dimensional  space  was  produced  directly  from  our  vortex-in-cell 
method,  (see  Chapter  IV,  section  2). 

Our  last  experiment  is  related  to  the  evolution  of  a  mixing  layer 
from  an  initial  state  of  uniform  vorticity  with  simple  two-  and  three- 
dimensional  small  perturbations.  We  expect  the  streamwise  perturbations, 
where  spanwise  uniformity  is  maintained,  to  lead  to  the  usual  roll-up 
of  vortex  patterns^^’^’^^ .  Combined  with  a  spanwise  perturbation  of 
the  same  type,  streamwise  distortions  of  the  vortex  filaments  occur. 
Again,  a  movie  tracing  the  vortex  filaments  in  three-dimensional  space 
was  produced. 

An  outline  of  the  remaining  chapters  of  our  work  is  as  follows. 

The  fundamentals  of  vorticity  dynamics  relevant  to  our  technique  are 
discussed  in  Chapter  II  and  a  description  of  the  computational  method 
is  given  in  Chapter  III.  In  Chapter  IV,  the  results  of  computational 
experiments  are  presented  and  discussed  while  conclusions  and  sugges¬ 
tions  are  given  in  Chapter  V. 
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Chapter  II 


BASIC  PRINCIPLES 

We  assume  unbounded  incompressible  flow,  fully  periodic  in  each 
of  the  three  dimensions.  Starting  with  the  incompressible  Navier- 
Stokes  equations  where  we  assume  no  external  forces,  the  equation  of 
motion  of  the  flow  field  is  given  by: 

vP  +  w2I?  (2.1) 

Dij  du 

where  the  "material  acceleration"  is  defined:  —  s  ^  +  u  •  . 

Rather  than  solving  (2.l)  on  an  Eulerian  mesh,  we  want  to  trace  the 
motion  of  the  vortex  filaments  in  the  velocity  field  these  filaments 
create.  The  collection  of  vortex  filaments  in  each  periodic  box  forms 
the  vorticity  field,  ui  =  V  x  u  .  We  therefore  solve  the  incompressible 
vorticity  equation,  obtained  by  taking  the  curl  of  (2.l): 

^|  =  U  •  Vi?  +  (2.2) 

By  using  the  equation  of  continuity, 

v  •  u  =  0  (2.3) 

we  find  that  the  velocity  field  can  be  determined  kinematically  from 

v2 u  =  -  v  X  U)  (2.4) 

[28] 

In  the  following,  we  use  the  large-eddy  simulation  approach1-  where 
one  calculates  precisely  the  large-scale  motions  and  models  approxi¬ 
mately  the  effects  of  the  finer  subgrid  scales.  The  separation  of 
the  large  and  small  eddies  mathematically  can  be  achieved  by  filter¬ 
ing^-  If  f  is  some  flow  field  containing  all  the  scales,  we 


7 


define  the  large-scale  or  resolvable-scale  component  of  f  to  be  a 
convolution  of  f  with  a  filter  function  G(r)  , 

ftf)  mfff  Ctf-*')  *(*')  dr  (2.5) 

In  the  present  method,  the  computed  scales  of  motion  are  assumed  to 
be  essentially  inviscid.  Any  actual  viscous  effects  are  on  a  subgrid 
scale  and  are  incorp'orated  into  the  filtering  procedure  described  in 
more  detail  in  the  next  chapter. 

From  Helmholtz's  theorem,  the  motion  is  purely  kinematic  and  the 
vortex  filaments  follow  material  lines.  Also,  by  definition  of 
lU  ,  V  •  uj  =  0  ,  and  Kelvin's  theorem  states  that,  in  an  ideal  fluid, 
the  velocity  circulation  F  around  a  closed  "fluid"  contour  is 

dr 

constant  in  time,  that  is  =  0  where 


r  = 


d  1 


-If 


dt 


(2.6) 


Here  A  is  the  cross-section  area  of  the  filament.  In  particular,  the 
effective  vorticity  field  uu  in  each  periodic  box  is  obtained  from 
the  actual  vorticity  it  by  filtering: 

iu(?,t)  =^^G(r-r'  )  uu(r',t)  dr'  .  (2.7) 

The  unfiltered  vorticity  u)  is  generated  by  the  space  curve  r^(£,t) 
as  follows, 

3(?,t)  T±f  6(?-^(§,t))  ^  d§  (2.8) 


and  is  highly  singular.  Here  §  is  a  parameter  which  traces  each 
filament  along  its  length  at  any  instant  in  time.  The  summation  is 
over  individual  vortex  filaments.  The  evolution  of  each  space  curve 


is  determined  from  the  continuous  velocity  field  by 


ar^S.t) 

at 


-fff  C(V?')  u(r',t)  d?' 


(2.9) 


with  u  determined  from  (2.4). 

In  summary,  equations  (2.7),  (2.8),  (2.9)  and  (2.4)  describe  the 
physics  of  vortical  flow  which  must  be  discretized  for  solution  on 
the  computer. 


Chapter  III 

DESCRIPTION  OF  COMPUTATIONAL  METHOD 
For  a  variety  of  reasons,  the  vorticity  field  and  other  fields 
are  conveniently  expressed  in  Fourier  space,  just  as  in  the  more 
successful  numerical  attacks  on  the  turbulence  problem  by  solution  of 
the  Navier-Stokes  equation^0"-^.  The  main  reason  the  fields  compo¬ 
nents  are  recorded  in  spectral  form  is  that  calculus  (differentiation 
and  integration)  translates  to  algebra  (multiplication  and  division) 
in  the  spectral  domain.  Of  course,  this  implies  periodicity  in  all 
three  dimensions. 

Working  in  Fourier  space  also  provides  additional  benefit  from 
the  control  one  obtains  over  the  filtering  operation:  the  convolution 
integral  in  (2.5)  becomes  a  product  in  Fourier  space.  Translating 
in  Fourier  space  the  basic  principles  given  in  the  previous  chapter, 
we  can  draw  a  computational  chain  of  operations.  Starting  with  a 
given  set  of  vortex  lines  with  their  circulation  T  ,  we  want  to  find 
their  displacement: 


The  numbers  refer  to  the  equations  of  the  previous  chapter.  The 
different  parts  of  the  scheme  will  now  be  explained,  namely  the 
numerical  modeling  of  the  vortex  filaments,  the  interpolations  that 
take  place  in  order  to  use  the  fast  Fourier  transform,  the  shaping 
of  the  vortices  and  the  time-stepping  procedure. 
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III.l  Filament  Modeling 

In  our  model,  one  describes  each  vortex  filament  by  a  succession 
of  closely  spaced  markers.  Considering  a  single  vortex  In  (2.8),  we 

have  _  _ 

iu(r)  -  vj  6  (r-r  (§  )  )  ||  d§ 

at  an  Instant  t  ,  where  T  Is  given  by  (2.6).  Taking  the  Fourier 
transform,  we  get 

tu(ic)  mJfJ e"lk  r  6  (*“?(§))  ||  d|  dr  (3.1) 

If  we  now  discretize  r  into  piece-wise  linear  sections, 

*^J.J-1  *  +  (1-5)'j_i  »  1  , 

then 

i(S)  -£  r/1^-?  )  d5  (3.2) 

j-1  0 


where  m  is  the  total  number  of  markers  describing  the  filament  and 
r  =■  r*  .  Integrating  (3*2)  and  letting 


£ 


2 
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we  obtain 


u)(k) 


/  -*  -* 

\  |  -Ik'r. 

\e  J 

26j 


2e1 


-it'fcfr.+r.  , ) 

e  J  J 
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(3.3) 


in  r*  ,  /-*  -*  V 

-ik^fr.+r.  .)  . 

E  (*,-*.  .)  e  1  I’1 

j=l  J  J_1  e. 


The  labor  of  evaluating  the  trigonometric  functions  for  all  sections 
of  the  filament  in  (3.3)  to  obtain  one  Fourier  component  $(£)  would 
be  prohibitive  if  each  vortex  required  at  each  time  step  the  evalua¬ 
tion  of  all  the  Fourier  harmonics  it  .  All  this  might  turn  out  to  be 
as  expensive  as  the  primitive  process  of  summing  up  all  direct  inter- 
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actions  among  the  vortices  .  It  is  more  efficient  to  first  dis¬ 
tribute  vorticity  onto  the  grid  according  to  an  interpolation  process 
and  then  perform  an  FFT.  Considering  equation  (3«3).  it  amounts  to 
two  things.  First,  we  should  be  able  to  express  sin  ej/®j  as  a 


linear  combination  of  e 


*lc  ^ 

1  where  ?  is  a  function  of  r.  and 


r^  ^  .  Secondly,  the  resulting  trigonometric  exponential  functions 
should  be  approximated  onto  the  mesh  in  Fourier  space  so  that  the  FFT 

can  be  performed.  As  of  the  substitution  of  sin  e  /e  >  we  went 

r  35*\ 

back  to  Eq.  (3.2)  and  evaluated  it  by  gaussian  quadratureL  to 
obtain, 

A 


m 


3<e>  -  r£ 

J-l  2 


X  <  e 


-ik 


•*[(1+3“*^  +  (l-3"^)rj_1J+ 


+(l+3'*)r 


J-l 


Indeed  this  is  equivalent  to  the  approximation, 

-ie  3’*  i€  3-^ 
sln  ej  e  ^  +  e  ^  cos 

6j  "  2 

(See  Appendix  I  for  derivation.) 


(i) 
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for  arbi- 


What  we  want  now  is  to  replace  the  pure  harmonic  e 
trary  r  by  an  approxlmant  to  be  evaluated  on  the  discrete  epatial 
mesh. 

III. 2  Interpolation 

In  each  of  the  three  dimensions,  we  use  quadratic  spline  Inter- 

ikx  ikn 

polation  to  approximate  e  in  terms  of  e  .In  particular. 
Ilex 

e  Is  represented  in  the  general  interval  n-^<x£n+^  as 
the  superposition  of  three  parabolic  arcs  as  follows: 


ikx 


■♦08, 


2  n+1 


(x-n-t^)2  +  g  ($-(x-n)2) 


(3.4) 

+  2  Vl 

This  representation  ensures  continuity  of  function  and  slope  between 

adjacent  intervals.  The  spline  coefficients  are  usually  chosen 

sc  as  to  force  agreement  between  function  and  approxlmant  at  the 

centers  of  the  range  in  which  each  quadratic  polynomial  is  used.  In 

our  case,  we  shall  inject,  at  this  point,  only  the  information  that 

ik 

the  spline  coefficients  gR  must  share  the  phase  shifts,  e  ,  per 

interval,  with  the  function  to  be  interpolated.  Thus  we  take 

gn  =  S(k)  eikn  where  S(k)  is  independent  of  n  .  (3.4)  then  becomes: 

ikx  cA.'J1  (  jl\2  ik  (n+1 )  h  ,  .2\  ikn 

e  -»  S(k)  -  (x-n+^)  e  +l£  -  (x-n)  le 

(3.5) 

1  /  .1  s2  ik(n-l)] 

+  g  (n+^-x)  e  v  'I 

or 

elkx  _>  S(k)|l  -  sin  |  +  ix  sin  kj 
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Note  that  the  choice  4/(3+cos  k)  for  S(k)  would  force  agreement 
between  function  and  approximant  at  x  =  0  • 


We  observe  that  the  representative  substitution  for  a  general 
interval  (3*5)  means  one  pretabulates  the  spline  coefficients 

^  S(k)  f(k)  e^n  ,  not  functions  values,  and  one  applies  the  spline 
k 

weights  to  the  entries  in  this  table  for  evaluation  of  local  function 
values.  The  indicated  summation  over  k  is,  of  course,  done  as  a 
fast  Fourier  transform.  But  one  does  not  transform  the  original 
spectrum  f(k)  ,  one  transforms  S(k)  times  the  spectrum. 

The  adjustment  factors  S(k)  must  still  be  specified.  As  men¬ 
tioned,  in  conventional  interpolation,  one  pins  the  interpolant  to 
the  function  at  specified  points:  the  mid-points  in  each  interval 
for  quadratic  interpolation.  This  is  appropriate  when  one  has  no 
more  information  on  a  function  than  the  values  at  integer  abscissae. 

In  our  case,  we  have  more  information  available. 

Since  we  know  the  function  to  be  approximated  throughout  the 
interval  in  principle  from  its  Fourier  harmonics,  we  can  make  a  less 
biased  choice  of  S(k)  .  We  can  minimize  the  mean  square  error  over 
the  interval  by  choice  of  S(k)  .  Using  standard  procedures,  one  finds 
that  S(k)  P(x)  deviates  from  e  with  the  least  mean  square  error 
when 

S (k)  -/Vx)  e_lkX  ixjf  ^  p2(x>  dx  (3.6) 


and  the  mean  square  error  is  then 


1 


h 

P(x)  e' 


ikx 


P2(x)  dx 
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Here  P(x)  is  the  polynomial  in  brackets  in  (3*5).  With  our  new 
choice  of  S(k)  ,  we  can  define  the  second  order  spline  as  that 
second  order  polynomial  which  joins  smoothly  (up  to  the  first  deriva¬ 
tive  continuous)  with  the  corresponding  polynomial  in  the  adjacent 
intervals  and  which  departs  from  the  given  function  with  least  mean 
square  error. 

The  numerator  in  (3-6)  is  recognized  to  be  the  Fourier  transform 

of  the  second  order  spline  when  the  trignometric  functions  in  P  are 

expressed  in  terms  of  complex  exponentials.  The  second  order  spline 

is  known  to  be  the  second  convolution  of  the  rectangle  function  with 

itself.  The  rectangle  function  (=  1  for  Jx|  <  ^  ,  =  0  for 

|x|  >  %)  transforms  to  ^  sin  j^3^  .  Thus  the  numerator  in  the 

2  k 

expression  (3*6)  for  S(k)  is  (—  sin  — )  .  As  for  the  denominator 

D  in  (3*6)  one  finds: 

2 

_  16  +  13  cos  k  +  cos  k 


,  2  k  .  2  .4k 

=  1  -  sin  2  +  15  Sin  2 


The  formula  for  the  mean  square  error  now  becomes 


“(i sin  i)/D 


The  square  root  of  the  mean  square  error  is  plotted  versus  k  in 
Figure  1  .  One  finds  that  for  small  k  ,  the  rms  error  is  of  order 


kJ  ,  specifically  the  rms  error  is 


k3  ~  .0058  k3 


Figure  1  .  Square  root  of  the  mean  square  error  k  for  nearest  grid  point,  linear 
and  quadratic  interpolation  of  a  pure  harmonic. 


As  mentioned  in  the  previous  section,  the  purpose  of  interpola¬ 
tion  is,  in  fact,  to  distribute  vorticity  onto  a  grid.  Figure  2 
illustrates  a  one- dimensional  model  where  three  quadratic  spline 
distributions  of  vorticity  with  different  total  amplitudes  are  shown 
to  have  their  centers  located  at  positions  a,  b,  and  c.  It  also  shows 
that  three  nearest  grid  points  to  each  of  the  centers  will  share  the 
vorticity  distribution  according  to  the  spline  function  weighting  on 
them.  In  other  words,  this  is  a  particular  way  to  distribute  a 
finite-sized  vorticity  on  its  neighboring  grid  points.  In  our  example, 
grid  points  N-2  to  N  will  share  the  vorticity  at  "a"  with  N-l 
getting  most  weighting.  Similarly  points  N-l  to  N+2  will  share 
the  vorticities  "b"  and  "c".  It  is  the  array  of  this  vorticity  distri¬ 
bution  that  will  be  transformed  to  calculate  the  potential  array  and 
velocity  field  array.  The  quadratic  spline  weighting  is  superior  to 
the  zero-order  weighting  (NGP  model)  and  first-order  weighting  (CIC 
model)  in  the  sense  of  creating  less  field-noise  and  resulting  in 
smoother  simulation  functions^^’ This  is  an  obvious  conclusion 
since  vorticity  is  now  distributed  among  three  grid  points  instead  of 
one  or  two  as  in  the  other  models  and  the  interpolated  distribution  is 
quadratic  rather  than  a  piece-wise  step  function,  or  first-order 
linear  function,  with  discontinuous  derivatives.  There  is  also  a 
reduction  of  aliassing.  In  three  dimensions,  the  three  nearest  grid 
points  in  each  dimension  (27  in  all)  will  share  the  vorticity  distri¬ 
bution  according  to  the  spline  function  weighting  on  them. 

The  sensitivity  to  grids  being  introduced  for  interpolation  pur¬ 
poses  was  tested  by  Buneman^  and  Wang^-1^.  Buneman^^  showed 

that  substantial  distortion  of  the  local  stream  function  can  result 
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from  linear  interpolation  if  the  vortex  is  located  at  other  than  the 
symmetric  positions  in  the  grid.  By  using  least-mean-square  fitted 
quadratic  splines,  it  was  shown  that  the  stream  function  is  virtually 
insensitive  to  the  grids  such  that  the  same  undistorted  stream  func¬ 
tion  will  appear  no  matter  where  the  vortex  is  located.  Figure  3 
illustrates  that  at  four  different  subgrid  positions,  the  first  for 
the  vortex  center  directly  on  a  grid  point,  the  last  for  the  vortex 
center  in  the  middle  of  the  cell  and  the  two  cases  in  between  for 
typical  unsymmetrical  vortex  positions:  the  stream  lines  of  the 
vortex  are  perfectly  circular  and  unaffected  by  the  existence  of 
the  grid  points. 

III. 3  Shape  Factor.  Subgrid  "Viscosity" 

It  is  intuitively  obvious  that  low-|ic|  harmonics  are  interpo¬ 
lated  by  a  certain  tabulation  mesh  better  than  high-|lc|  harmonics. 
Aliassing  sets  a  limit  at  kmax  *  n/A  for  each  component  of  i? 

(A  3  mesh  spacing):  any  harmonic  with  a  k-component  higher  than  this 
will  be  misinterpreted  by  the  interpolator  as  a  corresponding  lower 
harmonic  with  all  k-c<~aiponents  lying  within  the  interval  (-  tt/A  , 
tt/A)  .  In  effect,  the  .nterpolator  will  add  a  mixture  of  overtones 
to  the  approximation  of  a  pure  harmonic.  A  reduction  of  aliassing 
was  already  achieved  by  going  to  a  higher  order  interpolation  than 
linear.  Another  way  to  suppress  al lasses  is  to  de-emphasize  harmonics 
in  the  range  near  k 

max 

The  effects  of  a  finite  cut-off  of  the  spectrum  on  the  physics 
to  be  computed  is  an  important  question  separate  from  the  question 
of  interpolation.  Two  aspects  of  the  cut-off  problem  are  worth 
emphasizing. 
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Firstly,  a  sharp  cut-off  in  Fourier  space  is  undesirable  (no 
matter  how  perfectly  each  harmonic  is  evaluated)  because  it  surrounds 
the  objects  that  interact  via  the  field  with  halos  in  real  space. 

The  halos  decay  only  weakly  with  distance,  like  A/r  .  If,  instead, 
the  spectrum  is  brought  to  zero  more  smoothly,  say  at  least  parabo- 
lically,  such  halos  become  attenuated  more  strongly.  A  bell-shaped 
cut-off  factor,  say  G^(|lt| )  ,  suggests  itself.  Consequently,  even 
if  the  object  of  interpolation  studies  is  to  push  up  the  maximum 
usable  ]  Tc|  ,  it  should  concentrate  on  performance  in  the  range  of  low 
and  intermediate  lie!  ;  |lc| -values  near  k  become  irrelevant. 

It  is  important  to  distinguish  between  the  factor  G  (|Ic|  )  intro¬ 
duced  into  the  spectrum  by  such  shaping  and  the  compensating  adjust- 
ment  factors  S  (k)  which  improves  interpolation  performance.  The 
square  is  quoted  here  because  both  the  vorticity  deposition  and  the 
local  effects  of  the  velocity  field  call  for  interpolations. 

Secondly,  any  cut-off  factor  <5^(|Ic]  )  implies  a  shape  for  the 
interacting  vortices  which  is  the  Fourier  transform  of  6(|lc]  )  . 

The  dual  appearance  of  the  shape  factors  is  emphasized  by  the  fact 
that  the  interpolation  errors,  likewise,  appear  twice,  as  mentioned 
above.  So  the  shape  enters  the  field  interaction  process  twice, 
namely  both  when  the  velocity  field  harmonics  are  excited  by  the 
local  sources  (vortices)  and  when  the  velocity  field  harmonics  react 
back  on  them.  Therefore,  any  such  factor  in  the  spectral  domain 
(introduced  primarily  for  the  purpose  of  fitting  the  field  harmonics 
into  a  finite  computer)  can  be  interpreted  physically  as  vorticity- 
spreading. 

There  are  good  reasons  for  introducing  shapes  of  the  interacting 
elements  even  when  no  interpolation  is  used  at  all  and  spectral  data 
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are  evaluated  precisely,  without  grids  and  tabulations.  When  a  pro¬ 
blem  is  solved  entirely  in  the  spectral  domain,  non-linear  terms, 
i.e.  products  of  transformed  variables,  manifest  themselves  as  con- 
volutions  .  Here  the  cascading  into  higher  harmonics  (for  which 
there  is  no  room  in  the  computer)  has  to  be  suppressed  artificially 
by  some  cut-off  or  shape-factor.  Secondly,  the  interacting  "fluid 
elements"  in  the  real  world  are  usually  much  more  numerous  than  those 
that  can  be  accomodated  in  a  computer  with  its  peripheral  storage. 

Each  element  in  the  computer  stands  in  for  a  swarm  of  real  elements. 

It  should  therefore  be  given  a  spread. 

In  our  simulation  so  far,  we  use  a  cubical  mesh  system,  that  is 
Ax  =  Ay  =  Az  =  1  .  Therefore,  we  note  that  the  shape  factor  G(r ) 

should  be  isotropic  in  space:  it  knows  no  coordinate  axes.  Hence 

G(|ii|  )  is  isotropic  in  ic-space,  being  a  function  only  of  |  ic |  . 

2 

The  compensating  adjustment  factors  S  (k)  ,  on  the  contrary,  must 
be  functions  of  the  components  of  ic  .  Since  one  has  an  a  priori 

freedom  of  choice  as  regards  G(r)  or  <S  ( |  lc |  )  ,  one  is  at  liberty  to 

tailor  G(r)  so  that  it  de-emphasizes  the  poorly  interpolated  (badly 
aliassed)  harmonics  near  ]k]  =  n  .  Nevertheless  a  gaussian  profile 
G(r)  seems  to  be  more  readily  accepted  as  realistic:  for  a  viscous 

[4o]  [4i] 

vortex  ring  with  small  cross-section,  Tung  and  Ting  and  Saffman 

found  that  the  distribution  of  vorticity  across  the  core  is  gaussian. 

It  has  also  the  advantage  of  being  well  confined  both  in  real  and 

Fourier  space.  However,  with  a  finite  band-width  available  for  <5 ( | Ic |  )  , 

one  can  only  approach  a  gaussian  profile  since  the  finite  k- range 

always  causes  some  small  residual  halos.  A  fair  compromise  is  achieved 

by  taking  a  cubic  spline  profile,  tending  to  zero  cubically  as  k 

max 
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Is  approached: 


1  -  \  |S|2(l  -  if)  for  1E|  <  §  , 

tt  ' 

_  l£Lj  for  I  *  ^  *  n 


This  choice  is  quite  unrelated  to  the  use  of  splines  for  interpola¬ 
tion.  The  cubic  spline  profile  was  used  for  the  quantity  <5^(|lc|  )  , 
not  for  G(|k| )  .  The  transform  into  real  space  is 


This  is  positive  everywhere  and  has  very  small 
sidebands  (or  aliasses):  the  first  satellite  is  two  thousandth  of 
the  center  peak.  Both  this  <5^(|lc|  )  and  the  corresponding  profile 

in  real  space  look  very  much  like  gaussians.  In  fact,  by  using  the 

r  o£T\ 

Central  Limit  theorem1  ,  one  sees  that 


c2(|S| ) 


exp^-6j|L), 


Figure  4,  which  implies 


6(1*1) 


~exp|" 


(3.7) 


transforming  back  into  real  space,  we  get 

G(r)  =  exp i 


(3.8) 


A  spread-out  vortex  with  gauss ian  profile  is  the  result  of  vis¬ 
cosity  having  acted,  for  a  time  proportional  to  the  effective  area  of 
spread^1^,  on  an  initially  ideal  filamentary  vortex.  However,  the 
area  of  the  spread,  i.e.  the  "age"  of  the  vortices,  is  fixed  as  shown 
by  (3.8).  Therefore,  all  vortices  represented  in  this  manner  have 
the  same  "middle-aged"  spread! 
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To  sunmarize  the  interpolation  and  shaping  operations,  we  start 
with  a  vorticity  field  having  singularities  on  each  filament,  Eq. 

(2.8) .  The  use  of  a  grid  with  spacing  A  and  Fourier  transform 
methods  eliminates  the  singularities  since  only  harmonic  numbers  up 
to  tt/A  are  recognized  and  the  vortices  are  thus  automatically 
broadened  by  what  is  approximately  equivalent  to  a  rectangular  filter 
in  Fourier  space.  In  addition,  further  shaping  of  the  vorticity  dis¬ 
tribution  is  achieved  by  explicit  application  of  another  filter  whose 
representation  in  Fourier  space  is  G(£)  .  The  effective  vorticity 
field  uu  is  then  represented  by  a  filtering  operation  on  singular 

id  as  shown  by  (2.7)  where  G  now  Indicates  the  transverse  profile 
of  the  vortex  filaments.  We  choose  <2(lc)  to  be  almost  gausslan 
(bell-shaped,  coming  to  zero  smoothly  before  the  spectrum  would  cut 
off  more  drastically),  (Figure  4  )  and  then  each  filament  acquires  a 
transverse  profile  in  real  space  that  is  nearly  gausslan.  The  shape 
factor  G  is  applied  again  when  one  calculates  the  evolution  of  each 
space  curve  in  time  from  the  continuous  velocity  field  u(r,t),  Eq. 

(2.9)  with  u  determined  from  (2.4).  In  practice,  this  means  apply¬ 
ing  the  square  of  the  transform  of  G  as  a  filter  in  Fourier  space. 

A 

The  use  of  the  filter  G  provides  damping  of  the  high  wavenumber 

components  of  the  field,  damping  that  would  otherwise  occur  through 

rp£n 

subgrid  scale  dissipation  or  by  viscous  dissipation1 
III. 4  Solving 

The  advance  of  vortices  from  one  time  step  to  the  next  in  our 

O 

code  requires  the  handling  of  (6  field  .components  *  32J  mesh  points  ■) 
196608  field  components  and  about  32000  "vortex-markers".  If  we  now 
consider  Eq.  (2.4)  in  Fourier  space,  the  velocity  field  at  the  loca¬ 
tion  of  a  "vortex-marker"  is  obtained  by  weighting  the  entries  in  the 


25 


table  of  spline  amplitudes  with  the  spline  weights.  The  latter  are 
deduced  from  the  relative  position  of  a  vortex  in  Its  interpolation 
cell.  The  spline  amplitudes  are  obtained  from  the  velocity  harmonics 
by  first  multiplying  with  a  factor  S(kx)  and  two  similar  factors 
which  have  and  kz  in  place  of  k^  ,  then  calling  a  three- 

dimensional  FFT  on  the  resulting  array. 

The  same  factors  appear  again  when  the  displacement  of  the  kC^ 
vorticity  harmonic  is  calculated.  The  interpolation  can  be  done  for 
the  sum  of  all  the  harmonics,  and  the  table  into  which  one  interpo¬ 
lates  are  then  the  FFT's  of  the  harmonics  of  vorticity,  modified  by 
the  factors  insuring  best  mean  square  fit  in  each  dimension. 

In  going  from  the  table  of  spline  amplitudes  for  vorticity  to 
the  table  of  spline  amplitudes  for  velocity  harmonics,  one  therefore 
has  not  only  to  perform  a  forward  and  backward  FFT,  with  Eq.  (2.4) 
in  Fourier  space  in  between,  but  one  must  also  introduce  the  squares 
of  the  spline  fitting  factors  indicated  above. 

Similarly,  we  mentioned  that  any  vorticity  shaping  factor  should 
be  introduced  both  when  the  local  velocity  field  action  on  the  distri¬ 
buted  vorticity  cloud  is  evaluated  and  when  its  excitation  of  the 
vorticity  harmonics  is  accumulated.  In  both  cases,  one  could  perform 
a  convolution  in  real  space,  but  it  is  much  quicker  to  replace  this 
by  a  multiplication  in  Fourier  space.  The  transform  of  the  shape 
factor  is  therefore  introduced  squared  along  with  the  above  mentioned 
spline  fitting  factors  in  the  course  of  solving  the  equations  for 
the  velocity  field  in  Fourier  space.  It  is  convenient  to  introduce 
the  (squared)  shape  factor  along  with  the  inverse  Poisson  operator 
l/|it|2  .  In  our  32^  code,  there  are  onlj  256  possible  different 
values  of  |ic|2  In  the  sphere  |ic|  £  kmax  (64  values  for  a  16^  code), 
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so  any  function  of  { It J  can  readily  be  pretabulated.  The  schematic 
description  given  below  summarizes  the  algorithm. 


r  dr/d?  dr/dt 


III. 5  Time  Advancing 

The  choice  of  a  time  differencing  method  is  dictated  by  the 
trade-off  between  the  increased  cost  per  step  (in  computation  time 
and/or  storage)  of  hlgh-accuracy  methods  and  the  large  time  step 
allowed  by  such  methods.  When  the  choice  is  made  on  this  basis,  we 
are  interested  in  knowing  only  that  the  method  is  stable  fpr  the 
time  step  chosen.  There  is  no  reason  to  choose  a  method  with  extra 
stability  if  the  cost  is  higher.  Two  commonly  used  methods,  which 
are  second-order  accurate  and  require  only  one  function  evaluation 
per  time  step,  are  the  leapfrog  and  the  Adams-Bashforth  methods. 

Both  are  multi-step  explicit  methods  and  require  storage  for  two  time 
steps.  The  leapfrog  method  is  given  by 
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1 


n+1 

i 


+  2At  u" 


with  truncation  error 


and  the  Adams-Bashforth  method  is  given  by 

«  uuj  +  At(3u"  -  uJ_1)/2 

5  2  ^i 

with  truncation  error  -  ~  it  — ,  i  =  x,y,z  . 

12  btS 

f42  43] 

Both  the  methods  are  weakly  unstable1-  ’  .  Nevertheless,  used  with 

an  occasional  forward  Euler  step,  it  seems  to  be  possible  to  suppress 

T44  20I 

the  weak  instability  associated  with  leapfrog  differencing1-  ’  .  Also, 

a  linear  stability  analysis  shows  that  the  Adams-Bashforth  method 
exhibits  Improved  stability  over  the  leapfrog  method  and  that  its 

[44] 

total  spurious  computational  production  of  kinetic  energy  is  small 
In  this  work,  we  started  using  the  leapfrog  method  mainly  to  test  our 
code  on  the  rings  experiments  and  then  after  used  the  Adams-Bashforth 
method.  The  first  step  in  time  differencing  is  generated  by  Euler's 
method : 

=  U)°  +  At  u° 
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Chapter  IV 
RESULTS 

The  numerical  solution  of  the  scheme  described  in  the  previous 
chapters  for  different  initial  conditions  were  carried  out  on  the  CDC- 
7600  computer  at  NASA-Ames  Research  Center.  The  computing  time  per 
computational  time  step  to  move  a  vortex  made  of  m  markers  was  approxi¬ 
mately  as  follows: 

0.3^  +  m/5000  for  a  l6^  calculation, 

0.91  +  m/1275  for  a  32^  calculation. 

3 

In  our  32  algorithm,  some  programming  refinements  have  been  made;  since 
the  vorticity  field  and  velocity  field  are  purely  real  functions  in  real 
space,  it  is  possible  to  use  a  complex  Fourier  transform  of  length  l6  to 
get  the  vorticity  transform  of  length  32  and  similarly  to  transform  back 
the  velocity  harmonics  to  obtain  the  velocity  field  in  real  space.  (See 
Appendix  II  for  derivation  and^^.)  This  trick  improves  speed  as  well 
as  it  saves  storage. 

IV. 1  Single  Vortex  Rings 

A  first  experiment  was  done,  with  a  l6  mesh,  on  a  single  vortex 
ring  of  radius  R  about  the  z-axls.  Its  center  is  initially  located 
at  (8,8,8)  in  our  mesh  and  thereafter  moves  along  the  z-axis.  The 
circulation  is  I"  -  2  .  In  particular,  we  investigated  the  initial 
speed  of  the  vortex  ring  as  a  function  of  radius  and  position  around 
the  ring.  '» 

To  check  the  accuracy  of  our  mesh  technique  we  also  computed  the 
speed  of  the  vortex  ring  using  a  continuum  or  Green's  function  approach. 
Since  the  filter  we  use  in  the  mesh  method  is  approximately  gaussian 
(Eq.  3.8),  we  consider  a  single  vortex  ring  of  gaussian  cross-section. 
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On  the  one  hand,  if  we  consider  equations  (2.7)  and  (2.9)  in  Four¬ 
ier  space,  they  become  respectively. 


and 


ui(£,t)  =  G(|£|  )  •  uj(ie,t) 
iic*r . 


Sr  (5,t)  \ 3  rrr  ikt 

aT  ={h)Jffe  1  acl^l )  . 


(4.1) 


(4.2) 


On  the  other  hand,  equation  (2.4)  translates  to 


3(k,t)  -  lk-x 

k 


(4.3) 


in  Fourier  space.  Equation  (3.l)  together  with  (4.l)  gives,  for  a  single 
vortex  filament  of  gaussian  cross-section, 

ui(ic,t)  »  r  <5(1  ^1  )  J  e  lk  r  d£ 

Using  (4.2)  and  (4.3),  we  then  get 

r  (?,t)  -  — JSf  ii?  X  If’dl  rf 

at  (2V)3  J  k2  J 


(4.4) 


For  the  k-space  integration,  we  introduce  spherical  coordinates 

lc  *  |lc|  sin(9)cos(cp)e1  +  |£|  sin(9  )sin(cp)e2  +  |£|  cos^Je^ 
(see  Figure  5) 


Figure  5.  Vector  k  in 

spherical  coordinates. 


30 


31 


We  finally  obtain  for  (4.4), 

If  (5>t)  ■  fc. 1/0  gV-?) *  tf'd5  ^ 

where  2 

F(C)  =  erf  (C)  "  2tt_^C  e“C  (4.6) 

(4.5)  gives  the  filtered  velocity  of  translation  in  free  space  for  a 

single  vortex  filament  of  gaussian  cross-section.  It  should  be  mentioned 

2  2 

that  for  a  «  R  ,  (4.5)  can  be  approximated  to  yield 

f 

«  £ — flnf— Cl?  where  C  *  1.058  and  e  is  the  unit  vector  in 
3t  4ttRL  \ct  /  J  z  z 

the  direction  of  translation  z  (see  Appendix  III  for  derivation).  Note 
that  the  actual  speed  of  a  thin  vortex  ring  with  a  gaussian  distribution 
of  vorticity  has  been  calculated  by  Saffmarf^^  and  is  given  by  the  above 
formula  but  with  C  *  0.558.  The  difference  is  due  to  the  fact  that 
Saffman's  result  is  based  on  the  collective  motion  of  an  infinite  number 
of  vortex  tubes  with  internal  interaction  between  the  filaments  whereas 
our  result  represents  the  speed  of  a  single  computational  ring  filament. 

In  Figures  6  and  7 ,  we  plot  the  velocities  of  translation  measured 
at  two  positions  around  the  ring  versus  different  ring  radii  R  . 

Remember  that  the  vortex- in-cell  (VIC)  method  used  here  implies  periodic 
boundary  conditions  in  each  of  the  three  dimensions.  Therefore,  veloci¬ 
ties  measured  at,  say  (R+8,8,8)  and  (8+R/2^ ,8+R/2^ ,8)  are  not  quite  the 
same:  see  Figure  8.  In  order  to  explain  the  difference,  we  note  that 
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Figure  7  ,  Velocity  of  translation  maximum,  versus  radius  for 
periodic  array  of  single  vortex  rings. 


the  Fourier  transform  creates  images  in  a  cartesian  fashion  as  seen  in 
a  (x-y)-plane  (Figure  8).  The  images  oriented  diagonally  with  respect 
to  our  computational  ring  are  farther  so  they  will  not  slow  our  ring  down 
as  much  as  the  ones  on  the  axes.  Figure  6  shows  the  velocity  recorded  at 
points  of  the  ring  close  to  the  x-  and  y-axes,  where  the  velocity 
should  be  a  minimum  since  the  images  are  closer.  Figure  7  shows  the  velo¬ 
city  recorded  at  points  45  degrees  off  the  x-  and  y-axes,  where  the 
velocity  should  be  a  maximum. 

Our  results  are  compared  with  the  Green's  function  method  given  by 
(4.5)  to  which  were  added  the  Biot-Savart  contributions  of  the  periodic 
images.  (This  part  of  the  calculation  was  carried  out  using  Leonard's 
lagrangian  code  at  NASA-Ames^  Contributions  of  124  images  were 

supplied  in  this  manner.  Contributions  of  further  images  were  negligible). 
The  gaussian  width  used  in  the  Green's  function  calculations  was  chosen 

to  give  the  best  fit  to  the  vortex-in-cell  results  and  was  found  to  be 

2 

a  =  1.1  times  the  cell  area.  This  is  in  good  agreement  with  a  theore- 

2  2 

tical  estimate  of  cr  =  12/n  =■*  1.2  based  on  a  gaussian  fit  to  the  low 
| kj  behavior  of  our  filter  (Eq.  3.8).  Recall  that  our  filter  is  not 
strictly  gaussian  but  is  brought  smoothly  to  zero  at  |ic]  =  tt  . 

The  four  next  figures  show  pictures  of  the  initial  velocity  field 
in  the  middle  of  the  mesh  cells  for  a  ring  of  radius  R  =  4. 

Figures  9  and  10  show  the  field  in  the  planes  1.5  mesh  units  below 
and  above  the  plane  of  the  ring,  respectively  at  z  =  6.5  and  z  =  9«5  > 
where  the  magnitude  of  the  field  is  the  same  but  pointing  in  opposite 
directions,  respectively  toward  and  away  from  the  center  of  the  ring. 
Figures  11  and  12  show  the  field  in  the  planes  x  =  6.5  and  x  =  9*5  • 

Figure  13  shows  the  lateral  vortex  profile  in  the  (x-z)-plane  at  four 
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instants.  We  can  see  the  constancy  of  the  motion  in  time. 

A  second  test  was  done  on  a  set  of  two  vortex  rings  of  radius 
R  =  4  about  the  z-axis.  Their  centers  are  initially  located  at 
(8,8,7)  and  (8,8,10)  in  the  mesh  and  thereafter  move  also  along  the 
z-axis.  Both  have  the  same  circulation  r  «  2  . 

We  know  that  two  similar  vortex  rings  at  some  distance  apart  on  a 
common  axis  of  symmetry  will  do  the  following:  the  velocity  field 
associated  with  the  rear  vortex  ring  has  a  radially  outward  component 
at  the  position  of  the  front  ring  and  so  the  radius  of  the  front  ring 
gradually  increases  (with  T  constant).  This  leads  to  a  decrease  in 
its  velocity  of  translation,  and  there  is  a  corresponding  increase  in 
the  velocity  of  translation  of  the  rear  vortex  which  ultimately  passes 
through  the  larger  vortex  and  in  turn  becomes  the  front  vortex.  The 
maneuver  is  then  repeated.  Indeed,  we  observed  that  maneuver. 

Figures  l4, 15  and  l6  show  the  initial  velocity  field  respectively 
in  the  planes  z  *  5«5»  z  *  8.5  and  z  *  11.5.  As  expected,  at 
z  *  5.5  and  z  3  11.5,  the  magnitude  of  the  field  is  the  same  but 
pointing  in  opposite  directions,  respectively  toward  and  away  from  the 
center.  At  z  »  8.5,  centrally  between  the  two  rings,  the  field  reaches 
a  minimum.  Figures  17andl8  were  taken  at  x  -  6.5  and  x  =  9.5.  The 
last  five  figures  (I9-23  )  show  the  displacement  in  the  (x-z)-plane  and 
the  (x-y)-plane  at  ten  instants.  We  see  the  rings  going  through  each 
other  repeatedly  and  the  buildup  of  distortions  due  to  the  influence  of 


images 


6 


13 

10 

Y 

7 

4 

I 

I 

Figure  lU  . 


N  S, 


*  *  \ 

'  N  N 

"  N  \ 


>  >  t  i  t  t  ‘  ' 

\  \  l  J  l  J  ' 

\  \  11  l  /  '  '  ' 

W  \  \  t/S'  ' 

\\  \ / ///'  ' 


\ / / ^ 

. - -  ■>+%***’  « - - 


*  % 


/ 

t 

\ 

\ 

V 

/  X' 

/  / 

/ 

t 

\ 

\ 

\ 

V  ^  ^ 

✓  ^  s 

/  / 

/ 

t 

\ 

\ 

\ 

\  V  V 

*  *  / 

/  / 

/ 

1 

\ 

\ 

\ 

\  V  v 

•  *  * 

/  t 

t 

t 

1 

\ 

V 

X  %  v 

,  /  4  f 

*  lii  1 


*  \  \  \ 
—  »  -L  j — 


10 


13 


16 


Projection  on  (x-y)-plane  (z  =  5.5)  of  the  velocity  field  in 
the  middle  of  the  cells  for  a  pair  of  vortex  rings  (R  «  4) 
centered  at  (8,8,7)  and  (8,8,10)  respectively  in  a  l6xi6xi6  mesh. 
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Figure  16.  Projection  on  (x-y)-plane  z  ■  11.5)  of  the  velocity  field  in 
the  middle  of  the  cells  for  a  pair  of  vortex  rings  (R  ■  4) 
centered  at  (8,8,7)  and  (8,8,10)  respectively  in  a  l6xl6xl6  mesh. 
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Figure  17.  Projection  on  (y-zj-plane  (x  “6.5/  of  the  velocity  field  in 
the  middle  of  the  cells  for  a  pair  of  vortex  rings  (R  -  k) 
centered  at  3»S,7i  and  '8,8,10)  respectively  in  a  l6xl6xi6  mesh. 
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Figure  18.  Projection  on  (y-z)-plane  (x  =  9*5)  of  the  velocity  field  in 
the  middle  of  the  cells  for  a  pair  of  vortex  rings  (R  ■  U) 
centered  at  (8,8,7)  and  (8,8,10)  respectively  in  a  l6xi6xl6  mesh. 
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51 


IV. 2  Tavlor-Green  Vortex  System 

T241 

In  a  classic  paper,  Taylor  and  Green  considered  the  dynamical 
evolution  of  a  model  three-dimensional  vortex  field  in  order  to  clarify 
the  dynamics  of  turbulence.  The  relatively  simple  Taylor-Green  flow 
pattern  illustrates  the  basic  turbulence  decay  mechanisms:  production 
of  small  eddies  and  enhancement  of  dissipation.  The  latter  is  associa¬ 
ted  with  an  increase  of  vorticity  and  with  the  stretching  of  the  vortex 
lines.  The  Taylor-Green  system  has  also  proved  exceedingly  useful  for 
testing  numerical  and  perturbation  methods,  as  we  will  discuss  below. 

In  their  paper,  Taylor  and  Green  made  an  attempt  to  trace  the 
subsequent  motion  of  an  infinite  incompressible  fluid  with  velocity 
components  u^,  u^,  u^,  given  initially  by 

u^  =  A  cos (ax)  sin  (by)  sin(cz) 

Uy  =  B  sin(ax)  cos(by)  sin(cz)  (4.7) 

uz  =  C  sin(ax)  sin(by)  cos(cz) 

The  arbitrarily  assigned  constants  a,  b,  c.  A,  B,  C,  should  fulfill 
the  incompressibility  condition,  div(u)  =0  or  Aa  +  Bb  +  Cc  =  0  . 

The  flow  is  triply  periodic  and  a  typical  example  of  the  initial  flow 
pattern  for  a  given  set  of  constants  is  shown  in  Figures  24  and  25. 

At  an  early  stage  of  their  calculations,  Taylor  and  Green  realized  that 
it  was  impossible  to  obtain  significant  results  in  the  general  case 
when  the  initial  motion  is  represented  by  (4.7) »  so  they  confined  their 
attention  to  a  special  case,  namely  a-b*c,A*-B  and  C  =  0  . 

u^  ■  A  cos(ax)  sin(ay)  sin(az) 

Uy  *  -A  sin(ax)  cos(ay)  sin(az)  (4.8) 

u  -  0 
z 
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The  motion  remains  periodic  in  x,  y  and  z  for  all  time,  with  wave¬ 
lengths  2n/a  in  each  of  the  x,y  and  z  directions.  Taylor  and 

Green  used  power  series  in  the  time  t  to  find  u  ,  u  ,  u  ,  and  the 
r  x’  y*  z* 

pressure  p  at  time  t  .  They  took  averages  throughout  a  periodicity 

~ 2 

cube  and  deduced  the  series  for  the  mean  square  vorticlty,  uu  ,  as 
far  as  the  term  t'’  . 

In  order  that  the  effect  of  vortex  lines  stretching  may  be  illus¬ 
trated,  the  vortex  motion  must  be  three-dimensional  and  possess  a 
definite  scale.  It  is  difficult  to  think  of  a  motion  possessing  these 
characteristics  for  which  the  calculations  would  be  easier  than  for  the 
motion  (4.8).  That  is  the  reason  why  the  vortex  system  derived  from 

(4.8)  has  undergone  many  analytical  studies  since  the  original  paper 
[47-54  1 

of  Taylor  and  Green  .  Numerical  solution  of  the  Navier-Stokes 

equations  (2.1)  and  (2.3)  by  a  spectral  method  has  also  been  applied^0  ^ 
to  the  Taylor-Green  system  (4.8). 

In  our  case,  for  reasons  to  be  explained  below,  we  went  back  to 
(4.7)  and  considered  another  special  case,  namely  when  a  =  b  =  c  and 
-A  =  2B  =  2C  . 

u^  =  2A  cos (ax)  sin(ay)  sin(az) 

Uy  ■  -A  sin(ax)  cos(ay)  sin(az)  (4.9) 

u  *  -A  sin(ax)  sin(ay)  cos(az) 


See  Figures  24  and  25.  This  flow  is  symmetric  in  the  coordinates  y 
and  z  .  The  initial  vorticity  field  uu  =  V  x  u  is 


uux  -  0 


11)^  =  3«A  cos  (ax)  sin  (ay)  cos(az) 
uu^  *  3aA  cos  (ax)  cos  (ay)  sin(az) 


(4.10) 
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Projection  on  (y-z)-plane  (x  =  5*5)  of  the  velocity  field  in  the 
middle  of  the  cells  for  a  Taylor-Green  system  (eq.  4.9)  in  a 
I6xl6xl6  mesh. 


The  initial  vortex  lines  are  the  planar  curves 

sin(ay)  sin(az)  =  const, 
in  the  planes  x  =  const. 

See  Figure  26  for  an  illustration  of  the  vortex  lines. 

The  quantity  3  •  7u  ,  which  represents  amplification  and  rotation 

[igl 

of  the  vorticity  vector  by  the  strain  rate  ,  is  initially  nonzero. 
Therefore,  the  vortex  lines  given  by  Eq.  4.11  will  induce  a  velocity 
field  which,  while  transporting  the  lines,  may  also  change  their  lengths. 
Also,  since  the  initial  value  of  each  component  of  (u  X  J)  is  nonzero, 
the  three  components  of  velocity  continue  to  be  nonzero  after  the  initial 
instant  and  the  field  becomes  truly  three-dimensional.  (4.8)  and  (4.9) 
are  perhaps  the  simplest  examples  of  self-induced  vortex  stretching  in  a 
three-dimensional  velocity  field.  We  chose  to  study  the  system  (4.9) 
mainly  because  the  vortex  lines  (4.11)  (which  we  trace  numerically  in  our 
code)  are  easier  to  define  and  discretize  than  the  ones  resulting  from 
the  system  (4.8)  which  are  the  twisted  curves 

2 

sin(ax)/sin(ay )  =  const.,  sin  (ax)  cos(az)  *  const. 

On  the  other  hand,  we  do  not  know  of  any  paper  describing  a  perturba¬ 
tion  analysis  of  the  motion  initiated  by  (4.9).  Therefore,  before  going 
into  the  details  of  the  numerical  simulation  we  will  investigate  the 
evolution  of  our  Taylor-Green  vortex  system  (4.9)  by  developing  a  per¬ 
turbation  solution  to  the  Navier-Stokes  equations  (2.1),  (2.3)  in  powers 
of  the  Reynolds  number  Re  rather  than  in  powers  of  the  time  t  .  Gold- 

[47] 

stein  developed  such  a  series  solution  for  the  system  (4.8)  and  found 

that  this  series  is  more  inclusive  than  a  series  in  time  t  in  the  sense 
that  each  term  of  his  series  was  a  resummation  of  an  infinite  number  of 
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partial  terms  of  Taylor  and  Green's  series.  It  can  also  be  shown  that 
each  term  of  the  series  in  time  is  derived  from  only  a  finite  number  of 
terms  of  the  series  in  Re  . 


Hie  dynamical  problem  is  to  solve  the  Navier-Stokes  equations  for 
incompressible  flow: 

|^+u*  7u=-^Vp+  vV2!?  (2.1) 

v  •  u  *  0  (2.3) 

subject  to  the  (incompressible)  initial  conditions  (4.9).  The  pressure 
field  in  (2.l)  is  effectively  a  "Lagrange  multiplier"  that  insures 
compliance  with  the  incompressibility  constraint  (2.3).  Taking  the  diver¬ 
gence  of  the  equation  of  motion  (2.1)  and  applying  the  equation  of  con¬ 
tinuity  (2.3) ,  we  find  the  following  equation  for  the  pressure: 


v  •••  •&)’•  ef  &>' 


(du  du 
_J i 

ax  ay 


du  du  du  du 

+  _! l  —2 L  +  _ £ 
dy  dz  dz 


dx  / 


(4.12) 


We  reduce  all  the  equations  to  non-dimensional  form  by  writing 


u' 

=  u 

u' 

n 

c 

c 

II 

c 

T3 

'  =  £_ 

X 

X 

»  y 

»  2 

_Z  , 

? 

A 

A 

A 

pA^ 

x’ 

*  ax, 

y'  =  ay,  z'  = 

az,  t'  = 

Aat  , 

uo* 

=  uu 

Ul* 

=  U)  0)' 

=  u> 

X 

X 

.  y 

.  2 

2  9 

aA 

aA 

aA 

(4.13) 


Re  *  A 
va 
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From  then  on,  to  simplify  the  notation,  we  drop  the  prime  so  that  equa¬ 
tions  (2.3),  (2.1),  and  (4.12)  become 

v  •  u  =  0  (4.14) 


„  1  _2-» 
-  Vp  +  -  V  u 


•(&)'*&) 
+  2( 


Su  du  &u  &u 
_ Z _ x  +  _ z  _j 

dx  dy  Sy  3z 


+ 


The  Initial  conditions  (4.9)  are 


u  =  2  cos(x) 
x 

u  =  -  sin(x) 

y 

u  =  -  sin(x) 


sin(y)  sin(z) 
cos(y)  sin(z) 
sin(y)  cos(z) 


(4.15) 


(4.16) 


(4.17) 


at  t  -  0  . 

For  small  values  of  the  Reynolds  number  the  last  term  in  (4.15)  dominates 
while  the  quadratic  terms  remain  bounded  and  may  be  neglected;  thus  we 
have  to  solve  the  equation 


au 

at 


1 

=  —  ~u 
Re 


(4.18) 


with  the  initial  conditions  (4.17).  The  solution  is 


=  2  e"^t/,Re  cos(x)  sin(y)  sin(z) 

=  -  e"^t//Re  sin(x)  cos(y)  sin(z) 
yO 

A  *  -  e’^Re  sin(x)  sin(y)  cos(z) 

ZvJ 


(4.19) 
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For  the  second  approximation,  we  substitute  from  (4.19)  into  the 
quadratic  terms  in  (4.15)  and  into  (4.l6),  solve  (4.l6)  for  p  ,  sub¬ 
stitute  the  value  found  for  p  into  (4.15),  and  then  solve  (4.15)  again 
for  u  .  This  process  is  equivalent  to  substituting 


t/Re  =  T  , 


(4.20) 


expanding  u  , 


P 


in  series 

2 

-*  -»  -»  „  „  -> 
u  »  Uq  +  u^  Re  +  u^  Re  +  u^ 


3 

Re  + 


9 


2  3 

p  =  Pq  +  Re  +  p^  Re  +  Re  +  •  •  •  •  , 


(4.21) 


substituting  into  (4.15)  and  (4.l6)  and  equating  coefficients  of  powers 
of  Re  .  In  this  way  we  find  the  equations 


(4.22) 


+  2 


foiaijsi 

\ax  Sy 


au 

_s 

au 

du  _ 

au 

_2 

au 

rl  xO  .  zO 

i  +  — 

ax 

ay 

ay 

dz 

ay 

+ 


5u  _  3u  . 

—HQ.  — Sk  + 
az  ax 


au  au  n\ 
,  xl  _ z0» 

az  ax  ) 


and  so  on,  together  with 


6l 


(4.23) 


du 

^1  "  m  ?0  *  ^0  +  Vpo 


Su 


2->  p  ->  _»  _* 

^U2  •  ST”  =  u0  *  VU1  +  U1  *  %  +  Vpl 


and  so  on.  The  solutions  required  for  p  ,  p^,  etc.,  are  the  periodic 
solutions;  the  solutions  required  for  u^,  u^,  etc.,  are  periodic,  and 
are  zero  when  T  =  0  . 

Substituting  from  (4.19)  into  the  first  of  equations  (4.22)  we  find 

that 

V2p ^  =  \  e  ^4  cos(2x)  +  cos(2y)  +  cos(2z)  -  cos(2x)  cos(2y) 


-  cos(2x)  cos(2z)  -  4  cos(2y)  cos(2z) 


) 


(4.24) 


Hence 


PQ  *  -  e~^T  ^8  cos(2x)  +  2  cos(2y)  +  2  cos(2z)  -  cos(2x)  cos(2y) 

V 

-  cos(2x)  cos(2z)  -  4  cos (2y )  cos(2z)j  ^  ^5) 

From  the  first  of  equations  (4.23),  it  may  now  be  found  that 


V2uyi-  =  -  g  e'6T  ( cos  (2x )  sin(2y ) 


(4.26) 


Also  u  ,  =  0  when  T  =  0  .  Hence 

yi 


u  =  -  f^  cos(2x)  sin(2y) 


where 


.  3_  /  - 8t  -6t\ 

E1  =  16  \e  "  e  / 

_  Jl.lf  +  3Z_£  .  ITS*-  + 

"  3  \  8  8  +  12  24  +  • 


(4.27) 


(4.28) 
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uzl  ma^  ^oun<*  similarly.  It  is  easy  to  see  from  the  expressions 
for  Uq,  Pq  that  the  expression  for  may  be  obtained  from  that 


u^  by  interchanging  y  and  z 


Hence 

=  -  f^  cos(2x)  sin(2z)  (4.29) 

After  verifying  that  ux^  contains  no  term  independent  of  x  ,  we  find 
u^  from  the  equation  of  continuity,  the  result  being 

u^^  =  f^  (sin(2x)  cos(2z)  +  sin(2x)  cos(2y))  (4 . 30 ) 

Proceeding  in  the  same  way,  we  find  the  following  expressions  for 
p^,  u2,  p2,  u3  .  See  Tables  1,  2,  and  3-  We  note  that  the  expression 

for  u  „  is  obtained  from  that  for  u  by  interchanging  y  and  z  , 

z2  y  2 

and  so  on.  From  the  first  velocity  approximation  (i.e.  uQ)  on,  if  one 
expands  the  exponential  factors  (e’3\  the  g's,  h  )  in  Taylor  series, 
one  can  show  that  these  factors  are  ascending  powers  in  T  starting  res¬ 
pectively  from  the  zeroth  order  in  T  up.  Considering  (4.20)  and  (4.21 ), 
we  conclude  that  the  solution  of  the  equation  of  motion,  for  large  Rey¬ 
nolds  number  and  for  sufficiently  small  values  of  t  ,  can  be  represented 
by 


*+  -*(0)  / 
u  =  u^  '(x 

where 

un  (x,y,z) 

u  ,  n 
n 

-  0,  1,  2,  3 

expansion  of 
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TABLE  1 


Expressions  for  p  and 


sin(3x)  sin(y)  sin(z) 

sin(x)  sin(y)  sin(3z) 

sin(x)  sin(3y)  sin(z) 

cos(x)  sin(y)  sin(z) 

cos(3x)  sin(y)  sin(z) 

cos(x)  sin(3y)  sin(z) 

cos(x)  sin(y)  sin(3z) 

cos(3x)  sin(3y)  sin(z) 

cos(3x)  sin(y)  sin(3z) 

sin(x)  cos(y)  sin(z) 

sin(x)  cos(3y)  sin(z) 

sin(x)  cos(y)  sin(3z) 

sin(3x)  cos(y)  sin(z) 

sin(3x)  cos(3y)  sin(z) 

sin(x)  sin(y)  cos(z) 

sin(x)  sin(y)  cos(3z) 

sin(x)  sin(3y)  cos(z) 

sin(3x)  sin(y)  cos(z) 

sin(3x)  sin(y)  cos(3z) 


TABLE  2 


Expressions  for  p  and 


P2 

Ux3 

uy3 

cos (2x) 

m 

0 

0 

cos (2y ) 

l2 

0 

0 

cos (2z ) 

l2 

0 

0 

cos(2y)  cos(2z) 

h 

0 

0 

cos(2x)  cos (2y ) 

B9 

0 

0 

cos(2x)  cos(2z) 

Mm 

0 

0 

cos(2x)  cos(2y)  cos(2z) 

mm 

0 

0 

cos  (4x) 

x6 

0 

0 

cos (4y ) 

mm 

0 

0 

cos  (4z ) 

X7 

0 

0 

cos(4x)  cos(2y)  cos(2z) 

X8 

0 

0 

cos(4x)  co8(2y) 

l9 

0 

0 

cos(4x)  cos(2z) 

l9 

0 

0 

cos(4y)  cos(2z) 

110 

0 

0 

cos(2y)  cos(4z) 

110 

0 

0 

cos(2x)  cos(4y) 

X11 

0 

0 

cos(2x)  cos(4z) 

kl 

0 

0 

cos(2x)  cos(4y)  cos(2z) 

112 

0 

0 

cos(2x)  cos(2y)  cos(4z) 

l12 

0 

0 

sin(2x)  cos(2y) 

0 

-hl 

0 

sin(2x)  cos(2z) 

0 

-hl 

0 

cos(2x)  sin(2y) 

0 

0 

hl 

cos(2x)  sin(2z) 

0 

0 

0 

Uz3 

0 

0 

~0~ 

~ 

~0~ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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TABLE  3  --  Expressions  for  g  g  g  h  and  the 


We  see  that  when  is  expressed  as  a  triple  Fourier  series  the 

coefficients  of  the  various  terms  are  expressible  as  series  of  ascending 
powers  of  Re  multiplied  by  functions  of  T  ,  the  coefficient  of 
cos  x  8 in  y  sin  z  being  an  even  series  commencing  with  the  zeroth  power 
of  Re  ,  that  of  sin(2x)  cos(2z)  and  sin(2x)  cos(2y)  being  an  odd 
series  commencing  with  a  term  in  Re  ,  those  of 


cos(3x)  sin(y)  sin(z)  ,  cos(x)  sin(3y)  sin(z)  ,  cos(x)  sin(y)  sin(3z)  , 
cos(3x)  sin(3y)  sin(z)  ,  cos(3x)  sin(y)  sin(3z) 

2 

being  even  series  commencing  with  terms  in  Re  ,  and  so  on.  Similar 

results  hold  for  u  and  u  .  (See  equations  (4.19),  (4.27),  (4.29), 
y  z 

(4.30)  and  Tables  1-3.)  Hence,  going  back  to  the  "prime"  notation  and 
equations  (4.13),  if  we  write 


2  222  ,2  ,2  ,2  ,2  2,.2 

q  =  Ux  +  Uy  Uz  *  q  “  Ux  +  Uy  +  Uz  =  q  /A  * 


2  2  2  2  .2  .2  ,  ,2  .  ,2  2,  2.2 

U  =  (1)  +  U)  +  U)  ,  U)  =0)'  +  U)  +U)'  =  UJ  /a  A  , 

x  y  z  x  y  z 


(4.32) 


2  2  2 

and  denote  by  q'  and  uu'  the  space-average  values  of  q'  and 

p 

0) '  taken  throughout  a  periodicity  cube,  it  follows  from  the  orthogonal 

2  2 

property  of  the  terms  in  the  triple  Fourier  series  that  q*  and  u)’ 
will  be  expressible  as  even  power  series  in  Re  and  can  be  obtained 
correct  to  the  term  in  Re^  if  in  u^  the  coefficient  of 
cos(x')  sin(y')  sin(z' )  is  found  to  Re^  ,  that  of  sin(2x' )  cos(2z') 

and  sin(2x')  cos(2y')  to  Re  ,  and  the  remaining  coefficients  to 

P  2 

Re  ,  with  corresponding  results  for  u^  and  u^  .  Hence  to  find  q' 

2  4 

and  uj'  as  far  as  the  term  in  Re  it  remains  only  to  find  the 
coefficients  of  sin(x' )  cos(y')  sin(z')  in  u'y4  •  Further,  of  all 
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and 


-*  2 

the  terns  in  only  those  in  h^  are  required  in  finding  q' 

2  4  -* 

0) '  to  Re  :  the  remaining  terms  in  u^  (the  dotted  lines  in  Table  2) 

6  2  2 

affect  only  the  terms  in  Re°  in  q'  and  u)'  ,  and  need  not  to  be 

2  2  4 

found  if  we  are  finding  q'  and  uu'  only  up  to  Re 

In  the  same  way  when  p*  is  expressed  as  a  triple  Fourier  series, 
cos(2x')  ,  cos(2y')  ,  cos(2z')  ,  cos(2x')  cos(2y')  ,  cos(2x')  cos(2z')  , 
cos(2y')  cos(2z')  occur  with  coefficients  which  are  even  power  series 
in  Re  beginning  with  the  zeroth  power  of  Re  ,  sin(3x' )  sin(y' )  sin(z')  , 
sin(x')  sin(y')  sin(3z' )  ,  sin(x')  sin(3y')  sin(z' )  with  coefficients 
which  are  odd  power  series  in  Re  beginning  with  terms  in  Re  ,  and  so 
on,  as  can  be  seen  from  equation  (4.25)  and  Tables  1-3.  Hence  the 
average  values  of  the  pressure  and  of  the  pressure-gradient  are  expressi¬ 
ble  as  even  power  series  in  Re  ,  with  coefficients  which  are  functions 

of  T  ,  and  from  equation  (4.25)  and  Tables  1-3  the  constant  and  the 

2 

Re  term  can  be  found.  No  terms  in  sin(x')  sin(y')  sin(z’)  have  been 
found  to  occur  in  p^  .  No  further  calculations  of  the  expressions  for 
p’  have,  however,  been  carried  out,  so  the  average  values  of  the  pres¬ 
sure  and  the  pressure  gradient  could  be  found  only  as  far  as  the  term 
in  Re^  . 

The  coefficient  of  sin(x' )  cos(y')  sin(z' )  in  u^  was  found 
to  be  ,  where 


5315  e~3T  + 


3  e~9T  + 


1  141295616  16384 


-iZL 


-11T 


1376256 


991  e"15T  +  55163  e~17T 

450560  26492928 


+  401  T  e‘17T  -  44 1  e~19T  -  27  T  e"19T  -  3  e"2^T 

157696  2883584  180225"  3153920 


(4.33) 


The  expressions  for  u'  required  to  calculate  q'  and  u>*  to 
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the  terms  in  Re^  are  therefore  (see  equations  (4.19)»  (4.20),  (4.21), 
(4.27),  (4.29),  (4.30)  and  Tables  1-3) 

«  ^2e"3T  +  g1  Re2  -  2  1^  Re^j  cos(x')  sin(y')  sin(z') 

+  Re  -  hx  Re3)(cos(2y' )  +  cos  (2z' ))  sin(2x' ) 

+  11  ®2  1162  cos(3x')  sin(y' )  sin(z' ) 

"  11  82  R®2  (ain^y">  sin(z>  )  +  »in(y‘  )  sin(3z'  )jcos(x’  ) 

"  19  83  ^  (sin^3y'^  8in(z')  +  8ln(y')  sin(3z' )jcos(3x' )  , 

*  (-  e"3T  -  ~  gx  Re2  +  ^  Re^j  sin(x' )  cos(y’)  sin(z’) 

+  f1  Re  +  h^  Re3  j  cos(2x' )  sin(2y' ) 

"  11  82  R®2  sin(3x')  co8(y')  sin(z') 

+  11  g2  R®2  sln^x')  cos(3y')  sin(z') 

•  ^  g2  Re2  sin(x'  )  cos(y'  )  sin(3z’  ) 

+  Yg  83  Re2  sin(3x')  cos(3y' )  sin(z' ) 

uz  *  ("  e"3T "  ^3  8i 1162  +  ki  ***)  8ln^x'  )  8in(y'  )  cos(z' ) 

+  ^  Re  +  Re3  j  cos(2x')  sin(2z' ) 

-  §  gg  Re2  sin(3x' )  «in(y' )  cos(z'  ) 


69 


(4.34) 


O  p 

+  g2  Re  sin(x' )  8in(y')  cos(3z') 

”  11  g2  R*2  sin(x' )  sin(3y' )  cos(z' ) 

+  19  83  Re2  sin(3x')  sin(y' )  cos(3z’)  . 


Hence 


!(' 


-3T 


8,  ^ 


Re  - 


xiD4/4lS  2  ^  242  2  \ 

+  4  Re  \  n  82  +  36l83/  +  *  •  •  • 


+  Re 


F,  +  Re 


F2  + 


•) 


(4.35) 


where  F1  and  F2  are  given  in  Table  4. 

In  the  same  way,  if  we  write  down  the  formulae  for  tJ'  correspond¬ 
ing  to  (4.34)  and  take  the  mean  value  of  the  sum  of  the  squares,  we  find 
that 

U)'2  =  (e'3T  +  ^  g2  Re2  -  ^  Re4)  +  8  Re  -  h±  Re3) 

x  2  121  2  \ 

+  Re  \k  «2  +  18  83  )  +  *  ‘  *  * 

-  j*  (e'6T  +  Re2  Gx  +  Re4  G2  +  .  .  .  .)  (4.36) 


where  G^  and  G 2  are  given  in  Table  4. 

p 

There  is  some  independent  verification  of  these  results  since  q' 

2 

and  id'  are  related.  By  multiplying  the  equations  of  motion  in  order 

by  u  ,  u  ,  u  and  adding  we  obtain  the  usual  equation  of  energy,  which, 

x  y  z 

when  average  values  are  taken,  becomes 


(4.37) 
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TABLE  4 

Expressions  for  F^,  and  G 2 


F1 

F2 

G1 

G2 

-6t 

e 

-  1 

128 

.  6393 
70647808 

1 

'  128 

-6393 

70647808 

-12T 

e 

6? 

4096 

-1 

32 

E9 

H 

1 

0) 

-JL5 

128 

-  275 

344064 

--2 1 

128 

275 

147456 

-i6t 

e 

_3 

64 

59 

43008 

1 

8 

-59.— 

16128 

-l8T 

e 

0 

5741 

102400 

0 

17223 

102400 

-20T 

e 

0 

-  435971 
4139520 

0 

Te-20T 

0 

■  5867 

39424 

0 

■ggl 

-22T 

e 

0 

215465 
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(4.38) 

so  that 

(4.39) 

and 

a?  (*  ’2)  -  -  "“2  • 

(4 .40) 

i.e. 

d  u>'2 

dt'\%  q  /"  ■  r r 

(4.4l) 

Hence 

(4.42) 

From 

(4.35)  and  (4.36) 

it  follows  that  the  F  and  G 

are  connected 

by  the  relations 

Gl- 

1  **1  G  1^2 
’  6  dT  *  2  ’  6  dT  ’ 

(4.43) 

which  are  seen  from  Table  4  to  be  satisfied. 

Our  formulae  include  those  that  could  be  obtained  by  developing  in 
powers  of  the  time  t  .  The  trick  is  to  replace  the  exponentials  in  our 
formulae  by  their  expansions  in  T  and  then  use  (4.20)  to  get  the  per- 
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turbation  solutions  in  powers  of  t  .  Notice  that  finite  order  trunca¬ 
tions  of  (4.35)  or  (4.36)  do  not  have  the  singular  behavior  exhibited 


by  truncations  of  a  series  in  powers  of  t  as  t  -*  *°  .  However, 
examination  of  the  displayed  terms  of  (4.35)  or  (4.36)  suggests  that 
for  t/Re  ^  1  ,  perturbation  series  in  powers  of  Re  may  diverge  for 
Re  >  12  . 

Neither  perturbation  series  in  powers  of  t  nor  Re  can  describe 

the  evolution  of  the  flow  field  for  large  t  or  Re  .  In  addition  to 

the  fundamental  fluid  dynamical  interest  in  the  development  of  the 

Taylor-Green  vortex,  the  flow  is  a  most  convenient  one  on  which  to 

debug  and  perform  tests  of  our  vortex- in-cell  (VIC)  method.  For  these 

purposes,  we  compared  our  results  with  the  ones  obtained  from  a  program 

[55] 

developed  by  Rogallo  for  the  simulation  of  homogeneous  incompressible 

turbulence.  In  Rogallo' s  algorithm,  the  velocity  field  is  represented 

spatially  by  a  truncated  triple  Fourier  series  (spectral  method)  and 

followed  in  time  using  a  fourth-order  Runge-Kutta  scheme.  Notice  that 

for  the  Taylor-Green  vortex,  the  results  from  this  spectral  method  are 

infinite-order  accurate,  i.e.  errors  go  to  zero  faster  than  any  finite 

powe::  of  1/k  as  k  -»  00  ,  in  contrast  to  the  finite  order  accur- 
max  max  * 

[56] 

acy  of  difference  schemes  .  Since  there  are  no  experimental  results 
of  the  Taylor-Green  system  to  compare  with  and  the  perturbation  analysis 
being  woefully  inadequate  to  describe  either  the  large  t  or  the  large 
Re  behavior,  the  spectral  method  provides  the  best  results  to  compare 
with. 

The  Taylor-Green  system  is  one  of  continuous  vorticity  (eq.  4.10), 
and  it  is  important  to  represent  such  a  continuum  by  a  sufficient  number 
of  discrete  vortex  filaments.  These  are  given  by  eq.  4.11.  Assuming 
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our  computational  box  in  the  range  [-tt,  tt]  in  each  dimension.  Figure 
26  shows  the  initial  discretized  vortex  filaments.  Each  filament  having 
effectively  a  gaussian  core,  they  were  distributed  in  such  a  fashion  as 
to  fill  the  space  up  equally.  To  assign  the  proper  circulation  T  to 
each  of  them,  let's  consider  eq.  4.10  at  z'  =  ±  tt/2  (i.e.  z  =  ff/2a) 
which  implies  u£  =  (1)^  =  0  and  uT  =  -  3  cos(x')  cos(y')  .  Using  eq. 

2.6,  we  obtain 

r  .  -  zff™  (x1 )  cos  (y ' )  dx'  dy ' 

A 

where  the  surface  A  is  taken  to  be  a  square  centered  at  each  filament. 
This  gives  us  a  range  of  circulation  values  identical  for  all  octants  of 
our  space  except  for  their  sign  alternating  from  one  octant  to  another. 
Figure  24  and  25 display  that  feature. 

These  initial  conditions  being  set  in  our  code,  we  let  the  program 
run  and  made  up  a  movie  tracing  the  vortex  lines  in  the  three  dimensions. 
Figure  27  shows  four  stills  of  that  movie.  The  movie  was  taken  from  a 
l6xl6xi6  computational  box  with  the  vorticity  field  represented  by  l44 
vortex  filaments,  each  broken  into  an  average  of  24  nodes.  The  limita¬ 
tion  in  the  total  number  of  nodes  is  owing  to  the  3-D  graphic  system  that 
can  only  display  about  3700  three-dimensional  points  per  frame.  Also, 
the  graphics  package  draws  straight  lines  between  nodes  whereas  through¬ 
out  the  computations  sharp  corners  were  eliminated  by  smoothing  and  quad¬ 
ratic  spline  interpolation.  The  movie  showed  that  the  vortex  filaments 
themselves  do  not  remain  planar  and  that  there  is  a  considerable  tendency 
towards  disordering,  convoluting  and  stretching.  However,  it  was  observed 
that  the  symmetry  of  the  initial  conditions  in  the  y-  and  z-directions 
is  generally  preserved  throughout  the  entire  calculation  except  for 
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minor  discrepancies  which  occur  when  the  instabilities  make  the  flow 
erratic. 

Besides  being  able  to  follow  the  motions  of  the  vortex  lines,  our 
code  also  provides  automatically  all  kinds  of  spectral  quantities.  For 
later  purposes,  we  define  the  following: 


total  energy  =  ^*E(k)  dk  ■  %  q'2  =  ^  ^u^2  +  u^2  +  u^2j 
total  enstrophy  =J" k2  E(k)  dk  =  \  U)'2  =  ^  |u)^2  + 


streamwise  component  of  total  energy 


■/ E,<k> 


dk  =  \  u 


,2 


cross-flow  component  of  total  energy  =  ^E^(k)  dk  =  \  u^ 


■Aw 


,2 


C  2 

spanwise  component  of  total  energy  =J  E^(k)  dk  =»  %  u^ 


After  rescaling  our  problem  to  a  32X32X32  simulation  and  representing 
the  vorticity  field  by  1024  vortex  filaments,  each  broken  into  an  aver¬ 
age  of  30  nodes  (total  =  30720),  we  investigated  the  spectral  dynamics 
of  the  Taylor-Green  system. 

The  first  thing  we  did  is  to  run  a  movie  of  the  energy  spectrum. 
Figure  28  and  29  display  four  stills  of  that  movie.  On  a  semi-logari¬ 
thmic  scale,  the  energy  E(k)  contained  in  the  shells  of  radius  k  is 
plotted  versus  the  harmonic  radii  available  in  a  32X32x32  simulation. 

The  upper  right  corner  displays  the  running  time  and  the  total  enstrophy. 

The  seven  next  figures  (Figures  30  to  36)  display  the  energy  con¬ 
tained  in  the  seven  lowest  modes  as  time  goes.  Our  vortex- in-cell 
(VIC)  results  are  compared  with  the  ones  given  by  the  32x32x32  spectral 
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Figure  32.  Energy  in  the  third  harmonic  versus  time.  Vortex- in-Ce] 
method  compared  with  spectral  method.  32X32X32  mesh. 


8: 


Figure  3^.  Energy  in  the  fifth  harmonic  versus  time.  Vortex- in-Cell  (VIC) 
method  compared  with  spectral  method.  32x32x32  mesh. 


Figure  35.  Energy  in  the  sixth  harmonic  versus  time.  Vortex- in- Cell  (VIC) 
method  compared  with  spectral  method.  32X32X32  mesh. 


simulation  where  v  ■  0  .  Our  filter  d  was  chosen  as  to  fit  best  the 
spectral  results  at  t  ■  0  and  was  left  fixed  thereafter.  The  agree¬ 
ment  is  excellent,  specially  for  early  times. 

Figure  37  displays  from  top  to  bottom  the  total  energy,  the  stream- 
wise  component  of  total  energy  and  the  cross-flow  component  of  total 
energy  versus  time.  In  our  Taylor-Green  system, 

y*E  (k)  dk =/e  (k)  dk , 

y  Z 

by  symmetry,  as  we  pointed  out  earlier  in  the  perturbation  analysis. 

The  agreement  with  the  spectral  code  is  very  good.  Notice  a  slight 
increase  (~3*5$)  in  the  total  energy  at  t  =  7.2  .  This  is  an  indica¬ 
tion  that  the  stretching  and  the  fixed  core  filament  is  starting  to 
create  inaccuracy.  Recall  that  the  filter  d  is  fixed  from  t  =  0  on 
and  that  we  do  not  do  any  rediscretization  or  repacking  of  the  vortex 

filaments  after  the  beginning  of  the  calculation.  This  last  problem 

T  21  22  23] 

is  commonly  associated  with  Lagrangian  calculation1*  ’  ’  ,  that 

is  the  continual  addition  or  removal  of  particles  from  the  calculation. 
Depending  upon  the  physical  nature  of  the  problem,  particles  may  become 
crowded  and  yield  unrealistically  high  gradients  of  flow  variables  or 
the  number  of  particles  in  a  region  of  interest  may  be  so  low  that  no 
realistic  representation  of  the  flow  is  possible.  From  the  point  of 
view  of  maintaining  a  uniform  accuracy,  a  procedure  which  can  rearrange, 
add  or  delete  particles  as  necessary  should  be  applied. 

Finally,  Figure 38  shows,  on  a  semi-logarithmic  scale  the  total  en- 
strophy  versus  time  for  a  longer  time  span.  The  agreement  with  Rogallo's 
results  is  again  very  good  until  the  already  mentioned  problems  come  up. 
On  the  same  figure,  we  also  plotted  the  total  enstrophy  versus  time  when 
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Figure  37.  Total  energy,  x-component  and  y-component  of  total  energy 

versus  time.  Vortex-in-Ceil  (VIC)  method  compared  with  spectral 
method.  32X32X32  mesh. 
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Figure  38.  Enstrophy  versus  time.  Vortex-in-Cell  (VIC)  method  compared 
with  spectral  method.  32X32X32  mesh. 
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we  initially  set  an  arbitrarily  narrower  filter  G  for  the  entire  cal¬ 
culation.  This  was  like  broadening  the  vortex  profiles  or  making  the 
flow  more  viscous. 

In  conclusion,  it  is  clear  that  a  perturbation  analysis  of  the 
Taylor-Green  system  does  not  display  much  interesting  features  of  the 
flow  pattern  for  medium  or  large  t  or  Re  .  However,  our  vortex-in¬ 
cell  method  gives  very  good  agreement  with  a  purely  spectral  decomposi¬ 
tion  of  the  flow.  The  use  of  the  filter  (5  provides  damping  of  the 

high  wavenumber  components  of  the  vorticity  field,  damping  that  would 

[28] 

otherwise  occur  through  subgrid  scale  dissipation  or  by  viscous 
[ig] 

dissipation  .  The  nature  of  the  dissipation  provided  by  G  is  not 
well  established  at  this  time.  However,  based  on  comparisons  with 
spectral  calculations,  we  can  make  the  following  observations.  Invis- 
cid  simulations  of  three-dimensional  motions  which  have  no  scales 
smaller  than  the  grid  spacing  A  are  accurately  simulated  with  no 
energy  loss  by  the  present  method  using  a  broad  filter  <5  .  However, 
as  the  energy  cascades  to  scales  smaller  than  A  ,  some  dissipative 
effects  must  be  provided.  The  application  of  the  gaussian  filter  G 
provides  energy  dissipation  for  the  Taylor-Green  problem,  where  such  a 
cascade  exists.  As  a  final  observation,  we  pointed  out  that  a  reseeding 
procedure  should  be  incorporated  in  our  code  to  guarantee  numerical 
stability  at  later  times  in  the  calculations. 
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IV.  3  Shear  Laver 

In  many  flows  of  practical  interest,  there  are  interactions  between 
irrotational  regions  and  rotational  turbulent  regions.  As  an  example  of 
such  a  flow,  we  chose  to  simulate  the  mixing  layer.  In  a  mixing  layer, 
at  high  Reynolds  number,  the  regions  are  separated  by  a  thin  interface 
across  which  there  is  essentially  a  jump  in  the  vorticity  parallel  to 
the  layer  (see  Figure  39 »  showing  such  jumps  at  the  top  and  at  the  bottom 
of  the  layer).  The  usual  difficulty  in  simulating  such  a  flow  arises 
from  the  fact  that  such  a  vorticity  jump  would  be  diffused  by  finite 
difference  methods.  Using  the  vortex-in-cell  method,  we  can  retain 
sharp  gradients  in  vorticity  by  avoiding  difference  approximations  to 
the  term  representing  convective  transport  of  vorticity,  uij  •  Vu  (see 
Eq.  2.2). 

Large  coherent  eddies  have  been  observed  in  turbulent  shear 

layersL  ’  and  seem  to  play  an  important  role  in  the  mixing.  It  hat: 

[261 

been  observed1  that  the  pairing  of  large  eddies  is  central  to  the 
question  of  shear-layer  development.  Since  large  eddies  have  such  an 
importance,  it  is  of  considerable  interest  to  attempt  to  model  the  amal¬ 
gamation  process.  In  this  section,  we  follow  the  evolution  of  the 
mixing  layer  disturbed  by  simple  two-  and  three-dimensional  perturbations. 

In  our  coordinate  system,  as  defined  in  the  previous  section,  the 
x-direction  is  the  streamwise  direction,  the  y-direction  is  the  cross- 
flow  direction,  along  which  the  vorticity  jumps  are  encountered,  and  the 
z-direction  is  the  spanwise  direction,  along  which  the  unperturbed  flow 
is  uniform  (see  Figure  4o).  Our  mesh  system  is  32x32x32  cells.  As  said 
earlier,  using  fast  Fourier  transforms  to  solve  the  Poisson's  equation 
generates  periodic  boundary  conditions  in  each  of  the  three  directions. 
Therefore,  in  order  to  insure  continuity  of  the  velocity  field  in  the 
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cross-flow  direction  y  ,  we  cannot  attempt  to  study  a  flow  with  a 
single  velocity  shear  layer.  We  rather  investigate  a  parallel  flow 
with  two  opposite  linear  velocity  shear  layers,  that  is,  a  flow  where 
the  Initial  fluid  streaming  velocity  profile  is  an  even  function  with 
respect  to  the  y  axis  (see  Figure  4l).  The  initial  vorticity  field 
consists  of  two  uniform  layers  of  vorticity,  infinite  in  the  spanwise 
direction  2  and  having  opposite  vorticity  so  that  their  sum  equals 
zero.  In  our  simulation,  these  two  layers  are  replaced  by  arrays  of 
discrete  line  vortices  divided  evenly  into  two  groups  and  distributed 
separately  in  different  parts  in  space.  Those  with  positive  vorticity 
are  located  within  a  thin  layer  in  the  lower  half  plane  (from  y  =  -  9 
to  y  *  -  7)  while  the  others  carrying  negative  vorticity  are  distri¬ 
buted  in  the  layer  at  the  symmetric  position  in  the  upper  half  plane 
(from  y  =  7  to  y  *  9)-  That  is  to  say  y  ranges  from  “16  to 
+l6  (Figure  42).  This  distribution  ensures  the  maximum  spacing 
between  the  two  layers  and  their  periodic  images  in  the  cross-flow 
direction  y  .  The  results  shown  in  this  section  were  obtained  by 
initially  perturbing  only  one  of  the  layers  so  that  the  effect  of  the 
periodic  boundary  in  the  cross-flow  direction  y  can  be  monitored  by 
observing  subsequent  perturbations  by  the  initially  undistributed  layer. 
Each  layer  is  made  up  of  160  filaments,  and  each  filament  resolved  into 
32  nodes. 

As  a  first  case,  we  choose  a  streamwise  (x)  perturbation  of  the 
filaments:  y  «  yQ  +  sin(Trx/8),  x  ranging  from  -16  to+l6.  The  per¬ 
turbation  is  identical  for  all  z  (Figure  4o).  It  means  that  we  in¬ 
troduced  a  first-harmonic  down-stream  vorticity  variation  within  the 
layer  having  positive  vorticity.  This  sets  a  purely  two-dimensional 
problem  where  there  is  no  spanwise  z  dependence  for  all  times.  The 
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z-component  of  vorticity  is  the  only  non-zero  component^  °  .  Figure  40 
shows  the  initial  vortex  lines.  Figures  42  and  43  display  stills  of  the 
vorticity  contours  at  four  instants  in  time  and  taken  at  the  center  of 
the  computational  box,  z  =  0  . 

In  these  contour  maps  and  the  following  ones,  the  dashed  lines 
stand  for  negative  contours  and  the  solid  lines  for  positive  contours; 
the  tiny  numbers  along  the  contours  indicate  the  different  levels  and 
the  small  square  in  the  lower  right  corner  of  the  contours  indicate  the 
mesh  size.  It  should  be  pointed  out  that  the  graphics  package  draws 
straight  lines  within  the  mesh  size  which  explains  the  lack  of  smooth¬ 
ness  in  the  contours.  On  the  other  hand,  recall  that  our  algorithm 
which  governs  the  temporal  evolution  of  the  flow  guarantees  a  higher 
level  of  smoothness. 

Figures^  and  43  show  the  roll-up  into  two  vortices  as  has  been 
observed  by  other  investigators^1^’1  This  is  the  phenomenon 

of  amalgamation  referred  to  earlier.  Also,  the  initially  unperturbed 
top  layer  eventually  undergoes  the  two-dimensional  vortex  roll  distur¬ 
bances.  Notice  that  the  sense  of  the  pairing  rotation  in  the  positive 
layer  is  opposite  to  that  in  the  negative  layer.  No  pairing  of  these 
vortices  further  occurs  because  of  the  symmetry  of  the  initial  distur¬ 
bance.  Notice  also  that  during  the  initial  stages  of  the  rolling-up 
(up  to  TIME  ■  0.9) »  the  initial  structure  of  the  rolled-up  sheet  is 
still  evident.  However,  once  the  "eddies"  are  well  formed  and  begin 

to  rotate,  they  lose  this  structure  and  develop  into  isolated  blobs  of 

[14]  Ti41 

vorticity1-  .  Contrary  to  Acton's  simulation1,  ,  the  two  vortex 

blobs  at  TIME  «  3-6  will  not  merge  into  one,  due  to  the  symmetrical 

opposition  of  the  vortices  from  the  negative  layer  and  the  fact  that 

we  introduced  a  first-harmonic  variation  only. 
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Figures  44  and  45  display  stills  of  the  streamwise  component  of 
velocity  at  the  same  four  instants  in  time  and  same  position  in  space  as 
for  the  vorticlty  contours.  Similarly,  Figures^  and 47  show  four  stills 
of  the  cross-flow  component  of  velocity.  As  expected  in  this  two-dimen¬ 
sional  case,  the  spanwise  component  of  velocity  was  found  to  be  null 
(within  roundoff  errors).  Discussions  related  to  the  spectral  dynamics 
of  two-dimensional  perturbations  will  be  incorporated  with  the  ones  on 
the  three-dimensional  case  that  follows. 

Our  final  investigation  was  to  initialize  a  full  three-dimensional 
perturbation:  y  =  yQ  +  sin(nx/8)  •  sin(Trz/8),  x  and  z  ranging  from 
-16  to  +l6.  it  means  that  we  introduced  a  first-harmonic  streamwise  and 
spanwise  vorticity  variation  within  the  layer  having  positive  vorticity; 
the  filaments  as  well  as  the  nodes  making  each  filament,  both  followed 
a  sinusoidal  distribution  within  the  layer. 

This  problem  has  all  three  components  of  vorticity  and  velocity  pre¬ 
sent  for  t  >  0  .  Generation  of  streamwise  vorticity  is  observed  as  shown 
by  the  contour  plots  (Figures  50  and  51 )  and  the  streaks  formed  by  the 
vortex  lines  (Figures  48  and  49).  Figures  1^8  and  49  show  two  views  of  the 
vortex  filaments  of  the  perturbed  layer  evolving  in  time.  Here  we 
encountered  the  same  display  problem  as  with  the  Taylor-Green  system; 
owing  to  the  limitations  of  our  graphics  package,  only  one-third  of  the 
computed  nodes  could  be  shown  in  Figures  48  and  49.  That  is  the  reason 
why  we  refrain  from  displaying  all  the  vortex  lines  of  the  two  layers 
together  and  concentrate  attention  on  those  of  the  layer  initially  dis¬ 
turbed.  However,  in  all  contour  plots,  both  layers  are  displayed  to¬ 
gether. 

Quite  early  in  time,  the  display  of  the  vortex  filaments  becomes 


Figure  U5.  Continuation  of  Figure 
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obscure  due  to  the  intense  stretching  of  the  filaments  from  streamwise 
perturbations.  It  is  then  more  informative  to  look  at  the  different 


contours  of  vorticity  and  velocity.  The  twelve  next  figures  (50-61) 

show  six  groups  of  four  stills  for  different  contours  taken  at  the 

center  of  the  computational  box,  x  -  0  .  Sequentially,  we  have  the 

-z)  contours  of  u)  ,  U)  ,  (ti  ,  u  ,  u  ,  and  u  . 

x  y  z  x  y  z 

In  our  four  last  figures,  we  combine,  for  comparison  purposes,  the 
results  obtained  through  the  streamwise  perturbation  and  the  ones  ob¬ 
tained  through  the  three-dimensional  perturbation.  Figure  62  displays , 
on  a  semi-logarithmic  scale,  the  streamwise  and  spanwise  components  of 
total  energy  versus  time  for  both  cases.  Of  course,  the  two-dimensional 
mixing  layer  does  not  exhibit  a  spanwise  contribution  of  energy  since  its 
spanwise  component  of  velocity  is  null.  In  Figure  63  ,  we  see  the  cross- 
flow  component  of  total  energy  versus  time  for  both  cases;  since  the 


energy  variation  is  plotted  on  a  semi-logarithmic  scale  with  respect  to 


time,  the  initial  linear  section  of  the  curve  represents  the  exponential 
growth  of  the  component  energy  in  time,  and  its  slope  is  hence  the 
linear  growth  rate^®,'^’^,^‘^.  Figure  64  shows  the  total  energy  versus 


time  for  both  cases.  In  the  two-dimensional  case,  the  energy  stays  con¬ 


stant  whereas  the  three-dimensional  case  implies  stretching  of  the  vor¬ 


tices:  at  t  =  2.7,  the  energy  has  increased  by  only  one  percent  but 
at  t  =  3.6,  the  total  energy  is  up  by  eleven  percent.  Finally  in  the 
last  figure  (65),  the  total  enstrophy  versus  time  is  displayed  for  both 
cases. 


The  main  conclusion  to  be  drawn  from  these  computer  experiments  is 
that  the  third  dimension  is  extremely  important  in  the  evolution  of  a 
mixing  layer  even  though  it  is  an  "ignorable  coordinate"  in  the  idealized 
uniform  shear  flow. 
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MIN  =  -1 .0U9E+01  MOX  »  1.0U9G+01 


Y-VORTICITY 

TIME  =  .900  □ 


Figure  52.  Contours  uu  for  3-D  perturbation,  x  ■  0 
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Figure  53*  Continuation  of  Figure  52 
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Figure  5 U.  Contours  U)  for  3_n  perturbation,  x  *  0  . 
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Figure  56.  Contours  u^  for  3“D  perturbation,  x  ■  0  . 
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Figure  57.  Continuation  of  Figure  56. 
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Figure  64.  Total  energy  versus  time,  for  2-D( — )  and  3-D( - )  perturbations. 
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Figure  65.  Total  enstrophy  versus  time,  for  2-D(— )  and  3-D(— )  perturbations. 


Chapter  V 
CONCLUSION 

In  this  work,  we  have  developed  an  approach  to  three-dimensional, 
time-dependent  computations  of  flows  using  a  novel  vortex  tracing  method. 

More  than  thirty  thousand  finite-sized  vortices  were  advanced  in 
their  self-consistent  velocity  field,  using  a  second-order  explicit  time 
differencing  method.  The  vorticity  distribution  within  each  of  the 
vortices  is  given  a  gaussian  profile  by  means  of  a  shape  factor  imposed 
in  wave  vector  space.  A  machine-coded  Fourier  transform  is  used  to 
transform  into  and  out  of  this  space  at  each  time  step.  Poisson's 
equation  is  solved  by  dividing  by  the  square  of  the  magnitude  of  the 
wave  vector.  The  curl  operation  is  also  performed  in  wave  vector  space. 
Subgrid  accuracy  is  obtained  through  best-square-fit  quadratic  spline 
interpolations.  This  numerical  method  has  been  tested  and  no  undesir¬ 
able  grid  effects  or  numerical  instabilities  were  found. 

The  use  of  the  filter  6  provides  damping  of  the  high  wavenumber 
components  of  the  vorticity  field,  damping  that  would  otherwise  occur 
through  subgrid  scale  dissipation^®^  or  by  viscous  dissipation^"*"^. 

Based  on  comparisons  with  spectral  calculations  on  the  Taylor-Green 
problem  we  can  make  the  following  observations  concerning  the  nature  of 
the  dissipation  provided  by  (5  .  Inviscid  three-dimensional  flows  which 
have  no  scales  smaller  than  the  grid  spacing  A  are  accurately  simulated 
with  no  energy  loss  by  the  present  method  using  (5=1.  However,  as 
the  energy  cascades  to  scales  smaller  than  A  ,  some  dissipative  effects 
must  be  introduced.  The  application  of  a  gaussian  filter  (5  provides 
energy  dissipation  for  the  Taylor-Green  problem,  where  such  a  cascade 
exists,  but  produces  no  loss  of  energy  for  the  case  of  a  single  vortex 
ring  or  a  pair  of  coaxial  rings,  i.e.  when  there  is  no  cascade.  Based 


122 


on  this  limited  experience,  we  believe  the  use  of  a  filter  d  yields 
an  appropriate  model  of  subgrid  scale  dissipation. 

Compared  with  a  purely  Lagrangian  method  for  the  case  of  the  vortex 
rings  and  with  a  purely  spectral  method  for  the  Taylor-Green  system,  we 
saw  that  our  vortex-ln-cell  method  can  certainly  compete  in  computing 
time  efficiency  and  adequately  produce  the  same  accuracy  for  real  and 
spectral  dynamics  results.  It  is  a  great  advantage  to  be  able  to  gener¬ 
ate  direct  visualizations  of  the  real  flow  together  with  all  kinds  of 
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spectra,  so  useful  in  the  study  of  turbulence  .  This  advantage  was 
used  to  produce  movies  of  the  vortex  line  evolution,  for  the  Taylor-Green 
system  and  for  the  mixing  layer,  alongside  with  movies  showing  the  cas¬ 
cading  of  energy  from  the  lower  to  the  higher  modes. 

Also,  through  vortlcity  and  velocity  contours  taken  at  each  time 
step,  our  preliminary  study  of  the  mixing  layer  showed  the  now  well- 
established  two-dimensional  vortex  roll  disturbance.  The  exponential 
growth  in  the  cross-flow  component  of  energy  for  the  two-  and  three- 
dimensional  cases  was  also  observed.  It  was  remarkable  to  see,  in  the 
three-dimensional  perturbation  case,  distinct  streamwise  distortions 
and  "cusping"  due  to  the  stretching  of  the  vortex  lines. 

A  discussion  of  the  limitations  of  our  method  of  performing  fluid 
dynamical  experiments  on  the  computer  is  called  for,  namely: 

a)  the  fluid  is  treated  as  incompressible 

b)  gravitation  (buoyancy)  has  not  been  taken  into  account 

c)  the  viscosity  is  not  defined  explicitly 

d)  the  boundary  conditions  are  periodic 

e)  the  number  of  elementary  vortices,  although  large,  is  finite 

f)  all  elementary  vortices  are  of  the  same  finite-size  profile 


Limitations  (a)  and  (b)  can  be  overcome  by  including  consideration  of 
potential  flow  contributions  which  are  affected  by  pressure  and  buoy¬ 
ancy  (transport  of  vortices  is  not).  Our  method  could  then  be  applied 

T62  6^1 

to  studying,  for  instance,  Rayleigh-Taylor  instabilities1  *  . 

One  way  of  Introducing  explicitly  the  effect  of  viscosity  is  to 
change  the  size  of  the  elementary  vortices  according  to  their  age,  in 
other  words,  to  remove  the  limitation  (f)  .  Some  ad-hoc  methods  for 
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modifying  the  effective  size  of  vortices  have  been  suggested  by  Buneman1  J 
to  some  extent,  one  can  make  up  thinner,  younger  vortices  by  surrounding 
a  positive  (clockwise)  vortex  concentrically  with  a  halo  of  negative 
(anti-clockwise)  vortices.  However,  the  limit  of  resolution  is  set  by 
the  cut-off  imposed  in  ic-space,  and  such  narrowed  patterns  tend  to 
result  in  undesirable  sidebands  or  aliasses.  Chorin  and  Bernard^®^ 
treat  close  elementary  vortices  as  singular  filaments  by  direct  inter¬ 
action  and  let  only  distant  vortices  interact  via  the  grid.  By  close 
and  distant,  we  mean  less  or  greater  than  a  few  times  the  rms  vortex 
profile  radius.  They  then  use  a  random  walk  of  filamentary  vortices  as 
a  spreading  device  to  simulate  viscous  effects.  Their  work  was  two- 
dimensional.  Such  a  hybrid  method  (the  PPPM  method,  for  particle- 

particle/particle-mesh)  of  treating  interacting  objects  was  employed 

[641 

successfully  in  simulations  of  molecular  interactions1-  .  A  linked 
list  is  kept  which  allows  one  to  find  nearest  neighbors  up  to  any  chosen 
distance  and  calculate  their  direct  effect,  as  well  as  the  effect  which 
one  would  have  obtained  from  them  via  the  mesh.  The  difference  is  then 
added  to  the  total  mesh  field.  In  this  manner,  the  more  localized 
nature  of  the  flow  field  due  to  stretched  or  young  vortices  can  be  pro¬ 
perly  distinguished  from  the  broader  flow  field  of  compressed  and  old 
vortices. 
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Non-periodic  boundary  conditions  oan  be  simulated  in  a  variety  of 
ways.  One  way  is  to  use  readily  available  sine  and  cosine  transforms 
instead  of  complex  exponentials  to  simulate  (possibly  vortex- shedding) 
planar  walls  in  one  or  two  of  the  three  dimensions.  One  would  retain 
periodicity  in  the  third  dimension,  having  in  mind  the  simulation  of 
channel  flow  or  flow  through  a  re-entrant  wind  tunnel  of  rectangular 
cross-section  (but  with  curvature  effects  absent).  Another  variant  is 
the  simulation  of  "infinite"  boundary  conditions  using  similar  approaches 
as  in^^  and  in^^^.  In  the  case  of  the  mixing  layer,  the  use  of  per¬ 
iodic  boundary  conditions  is  justifiable  only  if  one  moves  with  the 
mean  speed  of  the  flow.  However,  the  size  of  the  eddies  grows  linearly 
with  the  streamwlse  distance,  and  one  reaches  a  point  at  which  the  size 
of  the  computational  box  must  be  increased.  This  problem  can  be  elimin¬ 
ated  if  a  way  can  be  found  to  work  in  an  infinite  domain.  From  a  prac¬ 
tical  point  of  view,  this  means  that  we  must  either  use  functions  appro¬ 
priate  to  such  a  domain,  e.g.,  Laguerre  or  Hermlte  polynomials,  or  use 
a  mapping  which  makes  the  transformed  domain  finite.  We  shall  look  into 
the  possibility  of  making  use  of  some  kind  of  mesh  mapping  in  the  cross- 
flow  direction.  In  order  to  preserve  the  applicability  and  usefulness 
of  the  Fourier  methods,  the  mapping  should  be  chosen  in  such  a  way  that 
these  methods  can  be  applied  with  little  loss  in  convenience  or  accuracy^* 
Restriction  (e)  is  merely  a  limitation  due  to  computer  memory  size 
and  will  be  alleviated  with  improved  computer  technology.  64x64x64  mesh 
codes  handling  almost  one  million  particles  are  already  becoming  avail¬ 
able  for  plasma  simulations  and  could  be  adapted  to  vortex  tracing. 
However,  the  large  data  bases  required  in  such  codes  call  for  non-tri- 
vlal  modifications  of  our  present  mode  of  operation:  first,  most  data 
will  have  to  reside  on  disc  most  of  the  time,  to  be  read  in  for  the 
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arithmetic  and  immediately  written  out  again.  Second,  a  faster  machine 
than  the  CDC  7600  is  needed  to  process  the  many  more  vortex  elements  in 
a  run  of  reasonable  duration:  one  thinks  of  a  Cray-1  or  the  soon  to 
come  Cray-2  computer! 

Eventually,  it  may  be  possible  to  treat  practical  flows  such  as 
airfoils,  combustion  chambers,  etc.,  by  our  method.  Before  that  can  be 
done,  much  more  effort  should  first  be  devoted  to  developing  an  explicit 
model  of  the  viscosity,  the  treatment  of  the  boundary  conditions  and 
the  mesh  layout  and/or  mapping. 


126 


Appendix  I 


NUMERICAL  INTEGRATION 


On  one  hand,  using  the  series  expansion  for  sine,  we  have, 

.  Crd\-i} 


letting  k  •  ~  -  , 


1  6  J  +  120  J  5040  j  + 


On  the  other  hand,  we  obtained  in  (3.2): 

m 


3(1?)  -  r£  (?  -?  )  i 

j-1  J  J 


/ 


where  1^  stands  for  /  e 
trapezoidal  rule  to  evaluate  Ij  ,  we  get: 


d§  .  By  using  the 


-ik'r  -ik*r 

ij  -  - - ¥■ - +  Error 


where 


Error 


hldid*  -ie,[??i+(1-5,?i-i] 


12 


0  <  §  <  1  3 


This  accuracy  to  the  first  order  in  «j  is  shown  in  the  resulting 
approximation  of 


-1.  i. 

-  e  J  +  e  J 
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*  COS  G j  *  1  “  2  ^  • 


To  increase  accuracy,  we  may  use  Simpson's  rule  to  evaluate  1^ 


/  -ilc'r  -ilc»%(r  +r  )  \ 

1^  =  i  le  J  +  4  e  J  J  1  +  e  J  1j+ 


Error 


where 


Error 


2880  o  <  §  <  1  “  l8o 


The  accuracy  is  now  to  the  third  order  in  as  shown  in 


sin  e. 


1  ,2 

3cos  el+3 


(i  .  !i_  +  iL  \  +  2 

\  2  +  24  •  •  •  •/  +  3 


_1  +  _L 

6  +  72  *  * 


We  have  increased  the  accuracy  of  our  approximation  but  we  needed  one 

more  exponential  so  it  increased  the  complexity  of  our  algorithm. 

[35] 

A  compromise  is  to  use  a  gauss ian  type  formula  for  Ij  , 
keeping  the  approximation  down  to  the  sum  of  two  exponentials.  Using 
the  zeros  of  the  second  order  Legendre  polynomial: 


P2(x)  =  J(3x2-1)  -  0 


which  implies  x 


i  ■  *ilf and  *2  ■  -VI  • 


we  obtain 
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Appendix  II 

PERFORMING  A  2N-POINT  TRANSFORM  BY  MEANS  OF  AN  N- POINT  TRANSFORM 
The  imaginary  part  of  a  complex  function  can  be  used  advantageously 
to  compute  the  transform  of  a  real  function  F(n)  defined  by  2N 


samples  by  using  a  discrete  transform  which  sums  only  over  N  values. 
That  is,  we  wish  to  break  the  2N  point  function  F(n)  into  two  N 
sample  functions.  We  divide  F(n)  as  follows: 


h(n)  =  F(2n) 
g(n)  *  F(2n+1) 


o=0,  1, 


N-l 


That  is,  function  h(n)  is  equal  to  the  even  numbered  samples  of  F(n), 


and  g(n)  is  equal  to  the  odd  numbered  samples.  (Note  that  h(n)  and 


g(n)  are  not  the  even  and  odd  function  decomposition  of  F(n).)  The 

[45] 

discrete  Fourier  transform  can  then  be  written  as 


#(k)  ® 


-i2nkn/2N 


-I;1  F(2n)  e-l«k(2»)/2»  +£  F(2wl)  e-l^k(2^1)/2» 
n=0  n=0 

-Xf  F(2n)  e-l2"t"/N  +  e"1"*'”  ^  F(2«-l) 
n=0  n=0 

h(n)  e-i2TTkn/N  +  e"iTTk/N  £  g(n)  e-12nkn/N 


n=0 

N-l 

22  8(n) 
n=0 


■i2nkn/N 


fi(k)  +  e"ilTk/N  g(k)  ,  k  =  0,  1,  ....  N-l 


(AII.l) 


To  efficiently  compute  fi(k)  and  g(k)  ,  let 


y (n)  =  h (n )  +  i  g(n)  . 

That  is,  y(n)  is  constructed  to  be  the  sum  of  two  real  functions  where 

one  of  these  real  functions  is  taken  to  be  imaginary.  From  the  linearity 

[45] 

property1"  ,  the  discrete  Fourier  transform  of  y(n)  is  given  by 
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y(k)  -  h(k)  +  i  £(k) 

-  (t^r(k)  +  i  fi^k))  +  i(gr(k)  +  i  g^k)) 

-  (*>r(k)  -  «t Ot) )  +  i(fit(k)  +  gr(k)) 

-  R(k)  +  i  l(k) 

where  we  used  the  decomposition  into  their  real  and  imaginary  parts  for 
li(k)  and  g(k)  .  We  decompose  both  R(k)  ,  the  real  part  of  y(k)  ,  and 
l(k)  ,  the  imaginary  part  of  y(k)  ,  into  even  and  odd  components'^ 

y(k)  =  %(R(k)  +  R(N-k))  +  %(R(k)  -  R(N-k) ) 

+  i  %(l(k)  +  l(N-k))  +  i  %(l(k)  -  I (N-k ) ) 

Knowing  that  a  real  even,  a  real  odd,  an  imaginary  even  or  an  imaginary 
odd  function  transforms  respectively  into  a  real  even,  an  imaginary  odd, 
an  imaginary  even  or  a  real  odd  function^^,  we  have 

h (k)  =  %(R(k)  +  R(N-k))  +  i  %(l(k)  -  i(N-k)) 

and  (All. 2) 

g(k)  =  %(l(k)  +  I (N-k) )  -  i  %(R(k)  -  R(N-k)) 

Substitution  of  (All. 2)  into  (AII.l)  yields 

F(k)  =  #r(k)  +  i  #.(k) 

where  the  real  part  of  the  discrete  Fourier  transform  of  the  2N  sample 
function  F(n)  is 

fr(k)  =  %(R(k)  +  R(N-k))+  Vrcos(nk/N)(l(k)  +  i(N-k)) 

-  l[sin(TTk/N)(R(k)  -  R(N-k))  (All. 3) 

and  similarly  the  imaginary  part  is 
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*iOO  *  *(l(k)  -  i(N-k))  -  %  sin(TTk/N)(l(k)  +  i(N-k)) 


-  %  cos (rrk/N ) (R(k)  -  R(N-k) )  (AII.4) 

These  formulas  hold  for  k  ■  0,  1,  N-l  but  since  F(n)  is  a  real 

function  defined  by  2N  samples  to  start  with.  Its  Fourier  transform  Is 

fog"] 

an  hermitlan  function1  J  and  the  rest  of  the  spectrum  Is  defined  from 
(All. 3)  and  (AII.4)  as  follows: 

fr(!H-k)  -  #r(N-k) 

k  =  0,  1 . .  N-l 

?i(N+k)  =  -  tfi(N-k) 


As  of  £(N)  ,  £^(N)  ■  0  since  F(n)  is  real  and  we  assign  f^(N)  to 
be  zero  through  our  shape  factor. 

In  a  straightforward  manner,  (All. 3)  and  (AII.4)  can  be  inverted  to 

yield 

R(k)  =  %(*r(l0  +  Fr(N-k))  -  ^cos  (nk/N )  (^  £  (k )  +  F^N-k)) 

-  Vsln(rrk/N)(#r(k)  -  *r(N-k)) 
l(k)  =  i(fi(k)  -  fi(N-k))-  ^sin(nk/N)(#1(k)  +  ^(N-k)) 

+  ^cos(TTk/N)(#r(k)  -  #r(N-k)) 


Considering  the  graphical  representation  of  a  complex  fast  Fourier 
r45] 

transform  ,  the  scheme  developed  above  is  equivalent  to  half  a  level 
in  addition  to  the  log^N  levels  of  the  FFT  over  N  points. 
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VELOCITY  OF  A  SINGLE  VISCOUS  VORTEX  RING 


In  equation  (4.5),  we  consider  a  single  vortex  ring  of  radius  R 

P  P 

and  of  gaussian  cross-section  a  .  We  assume  that  cr  «  R  . 


Figure  66.  Vortex  ring 
of  radius  R 


From  Figure  66,  we  say  that  |rj  -  R  and  that  |r-r'| 
of  the  chord  subtending  the  angle  9  ,  therefore  | r-r 1 
By  definition  of  a  vector  product,  we  have  then 


where  e  is  the  unit  vector  in  the  direction  of  translation  z  .  Also 


for  an  arc  of  6  radians,  the  length  §  *  R0  so  (4.5)  becomes 


Letting  C,  -  sin(cp),  dC  -  cos(«p)  do,  so 


1 


With  C  ranging  from  0  to  1  ,  we  can  consider  a  binomial  series 

2 

expansion  for  (l-C  )  such  that 

^  T  .  T  a  “  a  '  •  •  •  • 


C(l-C2)%  ■  c  2  8  8 


Hence  (AIII.l)  is  rewritten 


I 4  If  i F  '  *)' ac] 


(AIII.2) 


2  2 

In  the  second  integral,  we  use  the  assumption  that  a  «  R  so  F  can 
be  approximated  to  1  (see  Eq.  4,6),  which  implies 


But 


( •  ‘)dC  ■  IS 10  h1"2*) 

In  2 


So, 


St 


&[fwM  «"**]*■ 


Letting  5  *  2^  RC  and  using  (4.6), 


£  ”  sfejjf*  "'ll  "f  (5)  -  j  e'?  ) d?  +  ln  2]  s 


After  integrating  by  parts  on  the  first  term  of  the  integrand. 
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S  -  &f-fvwy  -jf 

<~h  »  *] 


2*  R/cj  2 

21nfg)  e~^ 


d§ 


-  erfl 


Letting  £  ■  ?  ,  the  Integral  In  the  above  expression  becomes 


♦  C  2R2  /a2 


*/  « -  j/’ dC .  4/ 

*0  ("C)%  0  («C)%  oD2 


Mfi) 


2R2/a2 

Since  *  0  ,  it  is  safe  to  say  that  the  second  integral  of 


dC 


2.2 


the  right-hand  size  is  of  the  order  e"2R  , 


so 


*  f  f  C  dC  -  %  (-  21n(2)-  Y)  +  o(e‘2R2/<j2) 

0  (ttO*  '  ' 


where  Y  is  the  Euler's  constant. 

On  the  other  hand,  using  an  asymptotic  expansion  for  the  error  func¬ 
tion,  we  have 


-2R 2/o2  3  -2R 2/o2 

erf|“l«  1  -  22 - r—  +  ^  e 


■W)- 


R(2n)‘ 


R3(32tt)* 


Again,  it  is  safe  to  say 


that  erf(is)«  1  -  o^1^) 


2  2 

asymptotically  for  a  «  R  . 
Finally,  we  obtain 


If  *  life[ln(V‘5)+  *(2  ln^+  y)-  1  +  ln(2)+  o(e"2R  /ct  )Jez  (AIII*3) 
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